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We study an ordinary differential equation (ODE) arising as 
the many-server heavy-traffic fluid limit of a sequence of overloaded 
Markovian queueing models with two customer classes and two ser- 
vice pools. The system, known as the X model in the call-center liter- 
ature, operates under the fixed-queue-ratio- with-thresholds (FQR-T) 
control, which we proposed in a recent paper as a way for one ser- 
vice system to help another in face of an unanticipated overload. 
Each pool serves only its own class until a threshold is exceeded; 
then one-way sharing is activated with all customer-server assign- 
ments then driving the two queues toward a fixed ratio. For large 
systems, that fixed ratio is achieved approximately. The ODE de- 
scribes system performance during an overload. The control is driven 
by a queue-difference stochastic process, which operates in a faster 
time scale than the queueing processes themselves, thus achieving a 
time-dependent steady state instantaneously in the limit. As a result, 
for the ODE, the driving process is replaced by its long-run average 
behavior at each instant of time; i.e., the ODE involves a heavy-traffic 
averaging principle (AP). 



1. Introduction. We study an ordinary differential equation (ODE) 
that arises as the many-server heavy-traffic (MS-HT) fluid limit of a sequence 
of overloaded Markovian X queueing models under the fixed- queue-ratio- 
with-thresholds (FQR-T) control. The ODE is especially interesting, because 
it involves a heavy-traffic averaging principle (AP). 

The system consists of two large service pools that are designed to operate 
independently, but can help each other when one of the pools, or both, 
encounter an unexpected overload, manifested by an instantaneous shift 
in the arrival rates. We assume that the time that the arrival rates shift 
and the values of the new arrival rates are not known when the overload 
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occurs. We want the control to automatically detect the overload. The FQR- 
T control is designed to prevent sharing of customers (i.e., sending customers 
to be served at the other-class service pool) when sharing is not needed, and 
automatically activate sharing when the system becomes overloaded due to 
a sudden shift in the arrival rates. 

This paper is the third in a series of five papers. First, In [15] we initiated 
study of this overload-control problem and proposed the FQR-T control; 
see [15] for a discussion of related literature. We used a heuristic stationary 
fluid approximation to derive the optimal control when a convex holding 
cost is charged to the two queues during the overload incident. Within that 
framework, we showed with simulations that FQR-T outperforms the best 
fixed allocation of servers, even when the new arrival rates are known. The 
stationary point of the fluid model was derived using a heuristic flow-balance 
argument, which equates the rate of flow into the system to the rate of flow 
out of the system, when the system is in steady state. 

Second, in [16] we applied the heavy-traffic AP as an engineering prin- 
ciple in order to justify the ODE considered here to describe the transient 
fluid approximation of the X system under FQR-T after an overload has oc- 
curred. We observed that the FQR-T control is driven by a queue-difference 
stochastic process, which operates in a faster time scale than the queue- 
ing processes themselves, so that it should achieve a time-dependent steady 
state instantaneously in the MS-HT limit, i.e., as the scale (arrival rate and 
number of servers) increases; see §3.1. We argued heuristically that the ODE 
should arise as the limit of a properly-scaled sequence of overloaded X-model 
systems, provided that the driving process is replaced by its long-run aver- 
age behavior at each instant of time. We performed simulation to justify the 
approximations . 

The present paper and the next two provide mathematical support. The 
present paper establishes important properties of the ODE suggested in 
[16]. The fourth and fifth papers prove limits. In [17] we prove that the 
fluid approximation, as a deterministic function of time, arises as the MS- 
HT limit of a sequence of X models; i.e., we prove a functional weak law 
of large numbers (FWLLN). This FWLLN is based on the AP; see [4, 8] 
for previous examples. In [18] we prove the corresponding functional central 
limit theorem (FCLT) that describes the stochastic fluctuations about the 
deterministic fluid path. 

We prove convergence to the ODE in [17] by the standard two-step pro- 
cedure, described in Ethier and Kurtz [5]: (i) establishing tightness and (ii) 
uniquely characterizing the limit process. The tightness argument follows fa- 
miliar lines, but characterizing the limit process turns out to be challenging. 
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Indeed, characterizing the limit process depends on the results here. Thus, 
the present paper provides a crucial ingredient for the limits established in 
[17, 18]. 

The AP makes the ODE unconventional. The AP creates a singularity 
region, causing the ODE not to be continuous in its full state space. Hence, 
classical results of ODE theory, such as those establishing existence, unique- 
ness and stability of solutions, cannot be applied directly. Moreover, existing 
algorithms for numerically solving ODE's cannot be applied directly either, 
since the solution to the ODE requires that the time-dependent steady state 
of the fast-time-scale process (FTSP) be computed at each instant. Never- 
theless, we establish the existence of a unique solution to the ODE, show 
that there exists a unique stationary point; and show that the fluid process 
converges to its stationary point as time evolves. Moreover, we show that the 
convergence to stationarity is exponentially fast. The key is a careful analysis 
of the FTSP, which we represent as a quasi-birth- and- death (QBD) process. 
Finally, we provide a numerical algorithm for solving the ODE based on the 
matrix-geometric method [10]. 

Here is how the rest of this paper is organized: The next two 
sections provide background. In §2 we elaborate on the X queueing model 
and the FQR-T control; that primarily is a review of [15]. In §3 we provide 
a brief overview of the MS-HT scaling and a heuristic explanation of the AP. 
In §4 we introduce the ODE that we study in subsequent sections. In §5 we 
state out main result, establishing the existence of a unique solution. In §6 we 
establish properties of the FTSP, which depends on the state of the ODE, 
and whose steady-state distribution influences the evolution of the ODE. 
In §7 we define the state space of the ODE, and prove the main theorem 
about existence of a unique solution. In §8 we establish the existence of a 
unique stationary point and show that the fluid solution converges to that 
stationary point as time evolves. In §9 we prove that a solution converges to 
stationarity exponentially fast. In §10 we provide conditions for global state 
space collapse, i.e., for having the AP operate for all t > 0. In §11 we develop 
an algorithm to numerically solve the ODE (given an initial condition) , based 
on the theory developed in the previous sections. We conclude in §12 with 
one postponed proof. 

Additional material appears in an appendix, available from the authors' 
web pages. There we analyze the system with an underloaded initial state, 
and show that the approximating fluid models then lead to our main ODE 
in finite time. We elaborate on the algorithm and give two more numeri- 
cal examples. We also provide a few omitted proofs. Finally, we mention 
remaining open problems. 



4 



O. PERRY AND W. WHITT 



2. Preliminaries. This section reviews the highlights of [15], starting 
with a definition of the original X queueing model, for which the ODE serves 
as an approximation. 

2.1. The Original Queueing Model. The Markovian X model has two 
classes of customers, arriving according to independent Poisson processes 
with rates Ai and A2- There are two queues, one for each class, in which 
customers that are not routed to service immediately upon arrival wait to 
be served. Customers are served from each queue in order of arrival. Each 
class-i customer has limited patience, which is assumed to be exponentially 
distributed with rate 9i, i = 1,2. If a customer does not enter service before 
he runs out of patience, then he abandons the queue. The abandonment 
keep the system stable for all arrival and service rates. 

There are two service pools, with pool j having rrij homogenous servers (or 
agents) working in parallel. This X model was introduced to study two large 
systems that are designed to operate independently under normal loads, but 
can help each other in face of unanticipated overloads. We assume that all 
servers are cross-trained, so that they can serve both classes. The service 
times depend on both the customer class i and the server type j, and are 
exponentially distributed; the mean service time for each class-i customer 
by each pool-j agent is 1/^ij. All service times, abandonment times and ar- 
rival processes are assumed to be mutually independent. The FQR-T control 
described below assigns customers to servers. 

We assume that, at some unanticipated point of time, the arrival rates 
change, with at least one increasing. We further assume that the staffing 
cannot be changed (in the time scale under consideration) to respond to 
this unexpected change of arrival rates. Hence, the arrival processes change 
from Poisson with rates Ai and A2 to Poisson processes with unknown (but 
fixed) rates Ai and A2, where Aj < mi/fc^, i = 1, 2 (normal loading), but 
Ai > Hi,iTUi for at least one i (the unanticipated overload). Without loss of 
generality, we assume that pool 1 (and class- 1) is the overloaded (or more 
overloaded) pool. The fluid model (ODE) is an approximation for the system 
performance after the overload has occurred, so that we start with the new 
arrival rate pair (Ai,A2). 

2.2. The FQR-T Control for the Original Queueing Model. The FQR-T 
control is based on two positive thresholds, fei^ and £2,1, and the two queue- 
ratio parameters, an d ^2,1 • We define two queue-difference stochastic 
processes £1,2 (*) = Qi(t) - n,2<22(*) and D 2 ,i = r 2 ,iQ2(*) - Qi(t). As long 
as Di t 2(t) < k\ t 2 and -D 2 ,i(i) < &2,i we consider the system to be normally 
loaded (i.e., not overloaded) so that no sharing is allowed. Hence, in that 
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case, the two classes operate independently. Once one of these inequalities is 
violated, the system is considered to be overloaded, and sharing is initialized. 
For example, if Di^{t) > k\^, then class 1 is judged to be overloaded and 
service-pool 2 is allowed to start helping queue 1. As soon as the first class-1 
customer starts his service in pool 2, we drop the threshold A;i ; 2, but keep 
the other threshold k% \. Now, the sharing of customers is done as follows: 
If a type-2 server becomes available at time t, then it will take its next 
customer from the head of queue 1 if D±2(t) > 0. Otherwise, it will take 
its next customer from the head of queue 2. If at some time t after sharing 
has started queue 1 empties, or 1)2,1 (*) = ^2,1 then the threshold &i |2 is 
reinstated. The control works similarly if class 2 is overloaded, but with 
pool-1 servers helping queue 2, and with the threshold &2,i dropped once it 
is crossed. 

In addition, we impose the condition of one-way sharing: we allow sharing 
in only one direction at any time. Thus, in the example above, where sharing 
is done with pool 2 helping class 1, we do not later allow pool 1 to help class 
2 until there are no more pool-2 agents serving class-1 customers. Sharing 
is initiated with pool 1 helping class 2 when Z?2,i(i) > &2,i and there are no 
pool-2 agents serving class-1 customers. And similarly in the other direction. 

Let Qi(t) be the number of customers in the class-i queue at time t, and 
let Zij(t) be the number of class- i customers being served in pool j at time 
t, i,j = 1,2. Let qi(t) and Zij(t) be the fluid approximations of Qi(t) and 
Zij(t), respectively. With the assumptions on the X system and the FQR- 
T control, the six-dimensional stochastic process (Qi(t), Zij(t);i, j = 1,2) 
describing the overloaded system becomes a continuous-time Markov chain 
(CTMC) (with stationary transition rates). 

Once sharing is initialized, the control makes the overloaded X model op- 
erate as an overloaded N model, and keeps the two queues at approximately 
the target ratio, e.g., if queue 1 is being helped, then Q±(t) n,2<32(i)- If 
sharing is done in the opposite direction, then r2 t iQ2(t) ~ Qi(t) for all t > 0. 
That is substantiated by simulation experiments, some of which are reported 
in [15, 16]. 

In addition to the thresholds k± t 2 and ^2,1, discussed above, the model 
also includes shifting constants K12 and ac 2) i. The shifting constants may be 
introduced after the threshold is dropped, because it may be dictated by 
the optimal ratio function in [15]. Let q* and z* •, i,j = 1, 2 denote the fluid 
steady state values of qi(t) and Zij(t). (We will show that a unique steady 
state, or stationary point, exists for the fluid approximation in §8 below.) If 
the optimal relation between the steady state fluid queues is gf = r\ 2 ( ?2+ K i,2 
for some K\ t 2 £ M, where r\ 2 denotes the fluid optimal ratio, (assuming 
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that pool 2 needs to help class 1), as is the case when the holding cost is 
separable and quadratic with non-zero constant and linear terms, then we 
use the shifted FQR-T control. Shifted FQR-T centers about ki^ instead 
at about zero. For example, if class 1 is overloaded, then every server takes 
his new customer from the head of queue 1 if Z)jj(t) > «i 2- Otherwise, it 
takes the new customer from the head of its own class queue. We call that 
control shifted FQR-T since it keeps the two queues at a fixed ratio, but 
shifted by the constant k± i. We can think of FQR-T as the special case of 
shifted FQR-T with «i )2 = 0. 

The beauty of the control is that, for large-scale service systems, FQR-T 
and shifted FQR-T tend to achieve their purpose; i.e., they keep the two 
queues approximately in fixed relation. In the stochastic system this means 
that the two-dimensional vector (Qi(t), Q^if)) evolves approximately as a 
one-dimensional process. In the fluid model this approximation becomes 
exact; We no longer need to consider the three-dimensional process x(t) = 
(qi(t), qzft), zi^it)), since it is enough to consider z\^{t) together with only 
one of the queues. The other queue is determined by the first via the state- 
space collapse (SSC) equation qi(t) = ri^q^it) + K^j, depending on which 
way the sharing is performed. In [17] SSC is shown to hold asymptotically 
in the MS-HT limit. 

3. The Many-Server Heavy- Traffic Fluid Limit. In this section 
we briefly describe the convergence of the sequence of stochastic systems to 
the fluid limit, as established in [17]. Without loss of generality we assume 
that class 1 is overloaded, and receives help from service-pool 2. (Class 2 may 
also be overloaded, but less than class 1, so that pool 2 should be serving 
some class- 1 customers.) 

3.1. Many-Server Heavy- Traffic (MS-HT) Scaling. To develop the fluid 
limit in [17], we consider a sequence of X systems, indexed by n (denoted by 
superscript), with arrival rates and number of servers growing proportionally 
to n, i.e., 

- A n _ m n 

(3.1) A™ = — > Aj and mf = — - — > mi as n — > oo, 

n n 

with the service and abandonment rates held fixed. We then define the 
associated fluid-scaled stochastic processes 



(3.2) Q«(t) = *i!l and Z^(t) = -^-, i,j = 1,2, t > 0. 



n 



n 
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For each system n, there are threshold kf 2 and k^i, scaled so that 



(3.3) 




The first scaling by n is chosen to make the thresholds asymptotically neg- 
ligible in MS-HT fluid scaling, so they detect overloads immediately when 
they occur. The second scaling by y/n is chosen to make the thresholds 
asymptotically infinite in MS-HT diffusion scaling, so that asymptotically 
the thresholds will not be exceeded under normal loading. It is significant 
that MS-HT scaling shows that we should be able to simultaneously satisfy 
both conflicting objectives in large systems. 

There are also the shifting thresholds nfj, arising from consideration of 
separable quadratic cost functions; see §2.2, but we do not specify their scale. 
If sharing is taking place, then at some time it was activated by sending the 
first class- 1 customer to service pool 2. We thus need only consider an d 
the weighted-difference process Di 2 (t) = QT(t) ~ r i2 ( 22W- Note that if 
—> oo, then D™ 2 -> oo as n -> oo. Hence, we redefine the difference 
process. Let 



where k ee k±2 and r = r\ 2 - 

With the new definition in (3.4), we allow K n to be of any order less than 
or equal to 0(n); in particular, we assume that K n /n —> k for < k < oo. 
There are two principle cases: k = and n > 0. The first case produces FQR 
(after sharing has began); the second case produces shifted FQR (shifted by 
the constant K n ). 

With the new process D n in (3.4), we can apply the same FQR routing 
rule for both the FQR and shifted FQR cases: if D n (t) > 0, then every newly 
available agent (in either pool) takes his new customer from the head of the 
class- 1 queue. If D n {t) < 0, then every newly available agent takes his new 
customer from the head of his own queue. 

3.2. A Heuristic View of the AP. The AP is concerned with the system 
behavior when sharing is taking place; i.e., when some, but not all, of the 
pool 2 agents are serving class 1. That takes place when q\ = rq2 + K. In that 
situation, it can be shown that the queue-difference process D n in (3.4) is an 
order 0(1) process, without any spatial scaling, i.e., for each t, the sequence 
of unsealed random variables {D n (t) : n > 1} turns out to be stochastically 
bounded (or tight) in M. That implies that D n operates in a time scale that is 



(3.4) 



K n )-rQ^(t), t>0, 
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different from the other processes Qf and Z™ 2 , which are scaled by dividing 
by n in (3.2). With the MS-HT scaling in (3.1), in order for the two queues 
to change significantly (in a relative sense, which is captured by the scaling 
in (3.2)), there needs to be O(n) arrivals and departures from the queues. In 
contrast, the difference process D n can never go far from 0, because it has 
drift pointing towards from both above and below. Thus, the difference 
process oscillates more and more rapidly about as n increases. Thus, over 
short time intervals in which X n remains nearly unchanged for large n, the 
process D n moves rapidly in its state space, nearly achieving a local steady 
state. As n increases, the speed of the difference process increases, so that 
in the limit, it achieves a steady state instantaneously. That steady state is 
a local steady state, because it depends on x(t), the fluid limit x at time t. 

To formalize this separation of time scales, we define a family of time- 
expanded difference processes: for each n > 1 and t > 0, let 

(3.5) D?(s) = D n (t + s/n), s>0. 

Dividing s by n in (3.5) allows us to examine what is happening right after 
time t in the faster time scale. For each t, a different process D™ is defined. 
For every t > and s > 0, the time increment [t,t + s/n) becomes in- 
finitesimal in the limit. A main result in [17] (Theorem 5.3) is that, for each 
t > 0, 

(3.6) D r t l = {Dl l (s) : s > 0} => D t {s) = {D t (s) : s > 0} in D, 

as n — > oo, where the limit Dt = {Dt(s) : s > 0} is the FTSP. 

For each n, the control depends on whether or not D n (t) > 0. In turn, 
the limiting ODE depends on the corresponding steady-state probability of 
the FTSP, 

(3.7) ti,2(x(*)) = lim P(D t (s) > 0) 

s— >oo 

which depends on x because the distribution of {D t (s) : s > 0} depends on 
the value of x(t) £ R s . 

4. The ODE. We now specify the ODE, which is the main subject of 
this paper. We assume that class 1 is overloaded, even after receiving help 
from pool 2. Hence both pools are fully busy and some pool-2 agents are 
helping class 1, so that = mi, Z2,i(t) = and zi^it) + 22,2 (*) — m 2- 

As a consequence, we only need consider z\ i among these four variables. 

We introduce an ODE to describe the evolution of the system state, which 
here is the vector x(t) = (qi(t), q2(t), zipfy)). The associated state space is 
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§ = [0,oo) 2 x [0,m2]. In particular, we consider the autonomous ODE 
(4.1) ' 

x(t) = (qi(t),q 2 (t),z 1>2 (t)) = = y(qi(t),q 2 (t),Z h2 (t)), t > 0, 

where \& : [0,oo) 2 x [0,7712] — > R 3 can be displayed via 
(4.2) 

qi(t) = Ai - mi/Ji,! - ni t2 (x(t)) [zi t2 (t)fj.i,2 + (jn 2 - z lt2 (t))fx 2>2 ] - 0i<?i(i) 
q 2 (t) = X 2 - (1 - 7ri j2 (x(t))) [(7712 - zi,2(*))^2,2 + ^1,2(^)^1,2] - 2 q 2 {t) 
h, 2 {t) = ni >2 (x(t))(m 2 - zi i2 (t))/x 2i2 - (1 - Ki i2 (x{t)))z lt2 (t)n lt2l 

with 7Ti i2 : [0, oo) 2 x [0, 771,2] — >• [0, 1] defined by (3.7) when q\ — rq 2 = k, 
fti, 2 (x) = 1 when q\ — rq 2 > k and Tri t2 (x) = when q\ — rq 2 < k. We also 
consider the associated initial value problem (I VP) 

(4.3) x{t) = *(x(t)), x(0) = w 

for W(x) in (4.1) - (4.2). 

5. The Main Result. The state space § is a subset of M 3 with the 
boundary constraints: q± > 0, q 2 > and < Zi j2 (t) < m 2 . The differ- 
ential equation for Z\ 2 prevents its boundary states and 7712 from be- 
ing active, because ^1,2 (^) = K\ j2 {x{t))m 2 Li 2i2 > when z± i2 (t) = and 
^1,2 W = (1 — T^i, 2 {x{t))m 2 iJLi j2 < when Zi t2 (t) = m 2 . However, the queue- 
length constraints can alter the evolution. In general, we can have q%{t) < 
when qi(t) = 0, which we understand as leaving qiit) fixed at 0. However, 
we are primarily interested in overloaded cases, in which these boundaries 
are not reached. Then we can consider the ODE without constraints. 

Recall that the shifting constant satisfies k > 0. We consider the restricted 
state space S = [k, 00) x [0, 00) x [0, m 2 ]. We thus avoid the transient region 
in which q\ < rq 2 + k with q 2 = 0, where qi(t) > and q 2 (t) < 0, but q 2 
remains at while q\ increases to the shifting constant k. The restricted 
state space, with q\ > k is shown to be the space of the fluid limit of the 
system in [17]. We will also show in Theorem 5.1 below that the ODE cannot 
leave this restricted state space. 

It is convenient to specify the conditions on the model parameters in terms 
of the steady-state formulas for the queues in isolation. For that purpose, let 
qf be the length of fluid-queue i and let sf be the amount of spare service 
capacity in service-pool i, in steady state, when there is no sharing, 7 = 1,2. 
The quantities qf and sf are well known, since they are the steady state 
quantities of the fluid model for the Erlang-A model (M/M/mi + M) with 
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arrival-rate Aj, service-rate mi and abandonment-rate 9{] see Theorem 2.3 
in [21], especially equation (2.19), and §5.1 in [15]. In particular, 

(5.1) qt = {Xl ~ and ^L-A)', i = i j2) 

where (x) + = max{x,0}. It is easy to see that gfsf = 0, i = 1,2. We thus 
make the following assumption, which is assumed to hold henceforth. 

Assumption A. 

(I) The model parameters satisfy 9i(qf — k) > Mi2 s 2- 
(II) The initial conditions satisfy x(0) 6 § = [n, oo) x [0, oo) x [0, 7712]. 

We now explain these assumptions. Clearly, a sufficient condition for both 
pools to be overloaded is to have sf = s% = 0, i.e., to have no spare service 
capacity in either pool in their individual steady states. However, if sf > 0, 
both pools can still be overloaded, provided that enough class-1 fluid is 
processed in pool 2. To have the solution be eventually in S, we require 
that 6i(qf — k) > ^1,2^2- This condition ensures that service pool 2 is also 
full of fluid when sharing is taking place, i.e., z\p, (t) + ^2,2 (*) — W2 for all 
t > (assuming that pool 2 is full at time 0). To see why, note that when 
service-pool 2 has spare service capacity (sf > 0), sharing will be activated 
if qf > K, because q% = 0. Now, the maximum amount of class-1 fluid that 
pool 2 can process, while still processing all of the class-2 fluid (so that qi is 
kept at zero), is ii\^.s%- hence, ^1,2^2 ^ s ^ ne minimal amount of class-1 fluid 
that should flow to pool 2. On the other hand, 9iqf = Ai — fJ>iirni is equal 
to the "extra" class-1 fluid input, i.e., all the class-1 fluid that pool 1 cannot 
process. Some of this "extra" class-1 fluid might abandon (if q\ > 0). The 
minimal amount of class-1 fluid that abandons is 9\n (but n can be equal 
to zero). 

We thus require that all the class-1 fluid, that is not served in pool 1, 
minus the minimal amount of class-1 fluid that abandons, is larger than 
Mi, 2 S 2- With this requirement, pool 2 is assured to be full, assuming that 
it is initialized full. (If pool 2 is not initialized full, then it will fill up after 
some finite time period; see the appendix.) 

Remark 5.1. (class 1 need not be more overloaded than class 2) In this 
paper we are interested in analyzing the ODE (4.2) as given. Hence, in 
Assumption A we do not assume that class 1 is more overloaded than class 
2; i.e., we do not require that qf — K > rq%. This extra assumption is not 
needed for our results for the specified ODE. In contrast, this assumption is 



ODE VIA AN AVERAGING PRINCIPLE 



11 



needed in order to show that the ODE holds as the fluid limit, with class 1 
receiving help; see Assumption 1 in [17]. 

We exploit Assumption A to show that the boundaries of § in R 3 play no 
role. 

Theorem 5.1. x(t) e S for all t > 0. 

We give the proof in §8.4 after the necessary tools have been developed. 
Our main result establishes the existence of a unique solution. 

Theorem 5.2. (existence and uniqueness) For any wq € 8, there exists a 
unique function x : [0, oo) — > § such that, (i) for all t > 0, there exist 5(t) > 
such that x is right- differentiate at t, differ entiable on (t,t + S(t)) and 
satisfies the IVP (4.3) based on the ODE (4.1) over [t,t + 5(t)) with initial 
value x(t), and (ii) x is continuous and differ entiable almost everywhere. 

Theorem 5.2 has two parts: First, there is (i) establishing the local exis- 
tence and uniqueness of a conventional differentiable solution on each inter- 
val [t,t + 5(t)), for which it suffices to consider a single t, e.g., t = 0. Second, 
there is (ii) justifying an overall continuous solution. 

We prove Theorem 5.2 in the next two sections. The proof is tied to the 
characterization of m t 2 in (4.2) and (3.7), and thus the FTSP Dt. We need 
to determine conditions for the FTSP Dt to be positive recurrent, so that 
the AP holds, and then calculate its steady-state distribution in order to find 
7Ti2- Moreover, we need to establish topological properties of the function 
^1,2, such as continuity and differentiability. We discuss the FTSP Dt next. 

6. The Fast-Time-Scale Process. Recall that the FTSP D t is the 
limit of D^ without any scaling (see (3.6)), where D™ is the time-expanded 
difference process defined in (3.5) associated with the queue-difference stochas- 
tic process D n = (Q™ — n n ) — rQV; in (3.4). Since there is no scaling of 
space, the state space for the FTSP Dt is the countable lattice {±j ± kr : 
j, k £ Z} in R. To see this, first observe from (3.4) that D n has state space 
{±j ±kr — K n : j,k £ Z}. Next, because of the subtraction in (3.5), Df has 
state space {±j ±kr : j,k £ Z}. Finally, because of the convergence in (3.6), 
the FTSP Dt has this same state space. 

6.1. The Fast- Time- Scale CTMC. We fix a time t and assume that we 
are given the value x(t) = (qi(t), ^(t), Zipif))- in order to simplify the 
analysis we assume that r is rational. That clearly is without any practical 
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loss of generality. Specifically, we assume that r = j/k for some positive 
integers j and k without any common factors. We then multiply the process 
by k, so that all transitions can be expressed as ±j or ±k in the state space 
Z. In that case, the FTSP D t = {D t (s) : s > 0} becomes a CTMC. 

Let A+ (m, x(t)), \+\m, x(t)), fi+\rn, x(t)) and /j,+ \m, x(t)) be the tran- 
sition rates of the FTSP Dt for transitions of +j, +k, —j and —k, re- 
spectively, when Dt(s) = m > 0. Similarly, we define the transitions when 



D t (s) = m<0: \ { 1> (m, x(t)), X ( _'(m,x(t)), ^'(m,x(t)) and n ( _'(m, x(t)). 



These rates are the limits of the rates of Df as n — > oo with X n (t) => x(t). 
First, for D t (s) = m £ (— oo,0], the upward rates are 



corresponding, first, to a class-1 arrival and, second, to a departure from the 
class-2 queue, caused by a type-2 agent service completion (of either cus- 
tomer type) or by a class-2 customer abandonment. Similarly, the downward 
rates are 



(6.2) tfL } (m,x(t)) = ii 1 , 1 z 1 ,i(t) + 6iq 1 (t), ^>{m, x(t)) = A 2 , 



corresponding, first, to a departure from the class-1 customer queue, caused 
by a class-1 agent service completion or by a class-1 customer abandonment, 
and, second, to a class-2 arrival. 

Next, for Dt(s) = m G (0, oo), we have upward rates 



corresponding, first, to a class-1 arrival and, second, to a departure from 
the class-2 customer queue caused by a class-2 customer abandonment. The 
downward rates are 

(6.4) 

fl^\m,x(t)) = + H\,2Z\,2{t) + /U2,2(m 2 - £1,2 (t)) + 0191 (t), 

^\m,x(t)) = X 2 , 

corresponding, first, to a departure from the class-1 customer queue, caused 
by (i) a type-1 agent service completion, (ii) a type-2 agent service comple- 
tion (of either customer type), or (iii) by a class-1 customer abandonment 
and, second, to a class-2 arrival. 




A w (m,x(t)) = Ai, 

\V)(m, x(t)) = fJ,i, 2 zi,2(t) + ^2,2{m 2 - zi )2 {t)) + 6 2 q 2 (t) 



(6.3) 



AW,x(*)) = A 



\ { l\m,x(t))=9 2 q 2 (t) 
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6.2. Representing the FTSP Dt as a QBD. Further analysis is simpli- 
fied by exploiting matrix geometric methods, as in [10]. In particular, we 
represent the integer- valued CTMC D t = {D t (s) : s > 0} just constructed 
as a homogeneous continuous-time quasi-birth- and- death (QBD) process, as 
in Definition 1.3.1 and §6.4 of [10]. In passing, note that the special case 
r = 1 is especially tractable, because then the QBD process reduces to an 
ordinary birth- and- death (BD) process. 

To represent Dt as a QBD process, we must re-order the states appropri- 
ately. We order the states so that the infinitesimal generator matrix Q can 
be written in block-tridiagonal form, as in Definition 1.3.1 and (6.19) of [10] 
(imitating the shape of a generator matrix of a BD process). In particular, 
we write 



f B 


A 








A 2 


Ax 


A 








A 2 


A x 


A* 








A 2 


At 



\ '■ J 

where the four component submatrices B,Aq,A\ and A 2 are all 2m x 2m 
submatrices for m = max{j, k}. In particular, These 2m x 2m matrices 
B,Aq, A\ and A 2 in turn can be written in block-triangular form composed 
of four m x m submatrices, i.e., 

(6.6) B A t) and A,,(f 

for % = 0, 1, 2. (All matrices are also functions of x(t).) 

To achieve this representation, we need to re-order the states into levels. 
The main idea is to represent transitions of Dt above and below within 
common blocks. Let L(n) denote level n, n = 0, 1, 2, . . . We assign original 
states 4>(n) to positive integers n according to the mapping: 
(6.7) 

(j){2nm + i) = nm + i and (j){(2n + l)m + i) = —nm — i + 1, 1 < i < m. 

Then we order the states in levels as follows 

L(0) = {l,2,3,4,...m,0,-l,-2,...,-(m- 1)}, 

L(l) = {m + 1, m + 2, . . . , 2m, — m, — (m + 1), . . . , — (2m — 1)}, 

With this representation, the generator-matrix Q can be written in the form 
(6.5) above, where A\ groups all the transitions within a level, Aq groups 
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the transitions from level L(n) to level L(n+1) and A2 groups all transitions 
from level L(n) to level L(n — 1). Matrix B groups the transitions within 
the boundary level L(0), and is thus different than A±. 

To illustrate, consider an example with r = 0.8, so that we can choose 
j = 4 and k = 5, yielding m = 5. The states are ordered in levels as follows 

L(0) = {1,2,3,4,5,0,-1,-2,-3,-4}, 

1(1) = {6,7,8,9,10,-5,-6,-7,-8,-9}, 

1(2) = {11,12,13,14,15,-10,-11,-12,-13,-14}, ... 



Then the submatrices B^, B\, A[ and A i , which form the block matrices 
B and Ai, i = 0, 1, 2, have the form in (6.12) below, where 

(6.8) a + = A? + A?> + + and a_ = aL 5) + \^ + ^ + ^ . 

(We solve a full numerical example with these matrices in §11.3.) 

Henceforth, we refer to Dt interchangeably as the QBD or the FTSP. 

6.3. Positive Recurrence. We show that positive recurrence depends only 
on the constant drift rates in the two regions, as one would expect. Let 5+ 
and 5_ be the drift in the positive and negative region, respectively; i.e., let 

8+(x(t)) = j U^ixit)) - $(x{t)j) + k (X^(x(t)) - $\x(t))) 

(6.9) ; ( ; ; 

6-(x(t)) = j [X b \x(t)) - fJ,v\x(t))\ + k (A W (x(t)) - ^(xtf))] . 

Theorem 6.1. The QBD Dt is positive recurrent if and only if 

(6.10) S-(x(t)) > > 5+{x(t)). 

Proof. We employ the theory in §7 of [10], modified for the continuous- 
time QBD. We first construct the aggregate matrices A = A + = 
Aq + ^ + A2 and A~ = A$ + A~{ + A2 . We then observe that the aggregate 
matrix A is reducible, so we need to consider the component matrices A + 
and A~ , which both are irreducible CTMC infinitesimal generators in their 
own right. Let v + and v~ be the unique stationary probability vectors of 
A + and A~ , respectively, e.g., with v + A + = and u + \ = 1. The theory 
concludes that our QBD is positive recurrent if and only if 

(6.11) v+A+Kv+A+1 and v~ A^ 1< v~ A^ 1. 

In our application it is easy to see that both v + and u~ are the uniform 
probability vector, attaching probability 1/m to each of the m states, from 
which the conclusion follows directly. □ 
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The alternative cases are simplified by the following relation: 

, g <*-(>(*)) - <*+(>(*)) = U + ^)(Mi,2^i,2 + (m 2 - zi, 2 )) 

> (j + k)m 2 (m,2 A ^2,2) > 0. 

Hence there are only two cases in which the drift does not point inward: (i) 
5 + (x{t)) > and S-(x(t)) > 0, (ii) 8-(x{t)) < and 5 + (x(t)) < 0. In both 
cases the behavior is unambiguous: In case (i), clearly iri^ix^t)) = 1; in case 
(ii), clearly 7Ti j 2(x(t)) = 0. 

6.4. Computing tti^ ■ When the QBD is positive recurrent, the stationary 
vector of the QBD can be expressed as a = {a n : n > 0} = {a n j : n > 
0, 1 < j < m}, where a n = (a+, a~) for each n, with a+ and a~ both being 
1 x m vectors. Then the desired probability ni t 2 can be expressed as 

00 m 00 00 

(6.14) vri i2 = J^^a+j = = ^a n l+, 

n=0 j=l n=0 n=0 
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where 1 denotes a column vector with all entries 1 of the right dimension 
(here mx 1), while 1+ represents a 2m x 1 column vector, with m l's followed 
by m O's. 

By Theorem 6.4.1 and Lemma 6.4.3 of [10], the steady-state distribution 
has the matrix-geometric form 

(6.15) a n = a R n , 

where R is the 2m x 2m rate matrix, which is the minimal nonnegative 
solutions to the quadratic matrix equation Aq + RA\ + R 2 Ai = 0, and can 
be found efficiently by existing algorithms, as in [10] (see §11 below). Since 
the matrices Aq, A\ and A2 have the block-diagonal form in (6.6), so does 
R, with submatrices R + and R~ . 

Since the spectral radius of the rate matrix R is strictly less than 1 (Corol- 
lary 6.2.4 of [10]), the sum of powers of R is finite, yielding 



J2R n = (i-RY 



n=0 

Also, by Lemma 6.3.1 of [10], the boundary probability vector ao in (6.15) 
is the unique solution to the system 

(6.16) a (B + RA 2 ) = and al = a (I - R)' 1 ! = 1. 

Finally, given the above, and using (6.14), we see that the desired quantity 
7Ti 5 2 can be represented as 

(6.17) tti.2 = a (/-i?n 1 l+. 

For further analysis, it is convenient to have alternative representations 
for tti^{x). Let the vector 1 have the appropriate dimension in (6.19) below. 

Theorem 6.2. {alternative representations for vri^) Assume that 5 +(x) < 
< 5-(x), so that the QBD is positive recurrent at x. (a) For r = 1, 

/ s o~-(x) 

(6.18) 7ri, 2 (x) 



5~(x) — 5 + (x) 

(b) For rational r, we have the sub-block representation 

c4(x){I - R+{x)Y l l 



(6.19) 7Ti 2 (x) = — c 

V 7 ' V 7 (4{x){I -R+{x))- l l + a^{x){I - R~(x))- l l 

where we choose ao(x) to satisfy ao(B(x) + R(x)A2(x)) = 0, renormalize to 
oto(x)l = 1, which corresponds to multiplying the original ao( x ) by a con- 
stant, decompose ao(x) consistent with the blocks as ao(x) = (ajj"(x), (x)). 
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Proof, (a) When r = 1, the FTSP D t = D t {x) evolves as an M/M/l 
queue in each of the regions Dt(s) > and Dt < 0. Thus, we can look at the 
system at the successive times at which Dt transitions from state to state 
1, and then again from state 1 to state 0. That construction produces an 
alternating renewal process of occupation times in each region, where these 
occupation times are distributed as the busy periods of the corresponding 
M/M/l queues. Hence, %i t 2(x) can be expressed as 

/ s E\T+(x)} 
(6.20) 7T lj2 (x) - 



E[T+(x)]+E[T-(x)Y 

where T + (x) is the busy period of the M/M/l queue in the upper region, 
while T~(x) is the busy period of the M/M/l queue in the lower region. 
By the definition of A, these mean busy periods are finite in each region. In 
particular, 

(6.21) E[T±(x)] 



p±(x)(l-p±(x)) ^(x)-\±(x) |«5±(x)|' 

where p (x) = X ± (x)/p ± (x), A + (x) and p + (x) are the constant drift rates 
up (away from the boundary) and down (toward the boundary) in the upper 
region in (6.3) and (6.4), depending on state x, while X~(x) and p~{x) are 
the constant drift rates down (away from the boundary) and up (toward the 
boundary) in the lower region in (6.1) and (6.2); e.g., \~{x) = pH\x) + 
P-\x) with j = k = 1 from (6.2). 

(b) We first observe that we can reason as in the case r = 1, using a regen- 
erative argument. We can let the regeneration times be successive transitions 
from one specific QBD state in level with Dt < to a specific state in level 
1 where Dt > 0. The intervals between successive transitions will be i.i.d. 
random variables with finite mean. Hence, we can represent tti^x) just as 
in (6.20), but where now T + (x) is the total occupation time in the upper 
region with Dt(s) > during a regeneration cycle, while T~(x) is the total 
occupation time in the lower region with Dt(s) < during a regeneration 
cycle. Each of these occupation times can be broken up into first passage 
times. For example, T + {x) is the sum of first passage times from some state 
at level to some other state in level 1 where Dt(s) > 0. The regenerative 
cycle will end when the starting and ending states within levels and 1 
are the designated pair associated with the specified regeneration time. The 
successive pairs (i,j) of starting and ending states within the levels and 
1 evolve according to a positive-recurrent finite-state discrete-time Markov 
chain. 
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Paralleling that regenerative argument, we can work with the QBD matri- 
ces, as in (6.17), but now using an alternative representation. Since 1++1 = 
1, where all column vectors are 2m x 1, we can apply the second equation 
in (6.16) to write 

= « (J-i?)^ 1 l + 
71-1,2 a (/- J R)- 1 l+ + aoU-^)" 1 l-' 

Then we can choose ao to satisfy ao(B + RA 2 ) = 0, renormalize to aol = 1 
(which corresponds to multiplying the original ao by a constant), decompose 
ao consistent with the blocks, letting ao = {olq,olq ), to obtain (6.19). □ 

With the QBD representation, we can determine when the FTSP Dt 
is positive recurrent, for a given x(t), using (6.10), and then numerically 
calculate 7Ti2- That allows us to numerically solve the ODE (4.1) in §11. 
We will also use the representations (6.17), (6.18), (6.19) and other QBD 
properties to deduce topological properties of nip. 

7. Existence and Uniqueness of Solutions. This section is devoted 
to proving Theorem 5.2. For the local existence and uniqueness in Theo- 
rem 5.2 (i), we will show that the function in (4.2) is locally Lipschitz 
continuous in Theorem 7.1 below. That allows us to apply the classical 
Picard-Lindelof theorem to deduce the desired existence and uniqueness of 
solutions to the IVP (4.3); see Theorem 2.2 of Teschl [19] or Theorem 3.1 
in [9]. Afterwards, in §7.2 we establish the global properties in Theorem 5.2 
(ii). 

7.1. Properties of ^ . We divide the state space S = [«, 00) x [0, 00) x 
[0,m2] = {(<7i, <?2i 21,2)} of the ODE into three regions: 

(7.1) S b = {q 1 - rq 2 = «}, S + = {q\ - rq 2 > k}, §~ = {qt - rq 2 < «}, 

with § = § 6 U §+ U S~. The boundary subset E> b is a hyperplane in the state 
space §, and is therefore a closed subset. It is the subset of § in which SSC 
and the AP are taking place (in fluid scale). In § b the function tti^ can 
assume its full range of values, < ir\ i(x) < 1. 

The region S + above the boundary is an open subset of S. For all x £ S + , 
^1,2(2;) = 1- The region S~ below the boundary is also an open subset of S. 
For all x G S~, 7Ti 2(2:) = 0. It is important to keep in mind that, in order 
for S~ to be a proper subspace of S, both service pools must be constantly 
full (in the fluid limit). Thus, if x £ §~, then z\ 1 = tci\ and z\ % + ^2,2 = m 2 
(but q\ and q 2 are allowed to be equal to zero). 
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It is immediate that the function ^ in (4.2) is Lipschitz continuous on S + 
and S - , because iti^ix) = 1 when x G 8 + , and vri i 2(x) = when x G 8~, so 
that ^ is linear in each region. However, ^ is not linear on § 6 . To analyze 
\l/ on S b , we exploit properties of the QBD introduced in §6. We partition 
§ 6 into three subsets, depending on the drift rates in (6.9). Let A be the set 
of all x G S b for which the QBD is positive recurrent, as given in (6.10); i.e., 
let 

(7.2) A = {x G S 6 | <5_(s) > > S+(x)}. 
Let the other two subsets be 

(7.3) A+ = {x G § b | 5+(x) > 0} and A" = {x G S b \ S-(x) < 0}. 

By the relation (6.13), there are no other alternatives; i.e., S 6 = AuA + UA~. 
Observe that Tti y 2{ x ) = 1 in A + , while tti^(x) = in A - . 

From the continuity of the QBD drift-rates in (6.9), if follows that A is 
an open and connected subset of S 6 . Hence, A can be regarded as an open 
connected subset of since § 6 is homoeomorphic to R+ x [0, 771,2]. Just as 
for the open subsets S + and S— , if the initial value is in A, then the ODE 
will remain within A over some initial subinterval. In contrast, the situation 
is more complicated in A + and A~ . 

There is potential movement out of the region S 6 only from the sets A + 
and A - . To understand what can happen, let d(x{t)) = qi(t) — rq2(t) and 
d'(x(t)) = qi(t) — rq2(t), from (4.2). On A + and A - , the possibilities can be 
determined from the following lemma. 

Lemma 7.1. On S b , if 7ri j2 (x) = 1, then d'{x) = 5 + (x); if 7ri j2 (x) = 0, 
then d'(x) = 5_(x). Hence, on A + , d'(x) > 0, while on A~ , d'(x) < 0. 

Proof. Substitute the appropriate values of ^1^(^(7;)) into (4.2) and 
compute 5±(x) from (6.1)-(6.4), recalling that r = j/k. □ 

We next separate equality from strict inequality for the weak inequalities 
in Lemma 7.1. For that purpose, we decompose the sets A + and A~ by 
letting 

A+ ={xgA+ I 6+{x) > 0}, A+ = {x G A + I 5 + (x) = 0}, 

(7.4) AI ={xgA~ I 5_(x)<0}, Aq ={xeA- \ 5_(x) = 0}. 

Lemma 7.2. Suppose that a solution exists for the ODE over a suffi- 
ciently small interval starting at x(0). If x(0) G A^, then x(0+) G S + ; if 
x(0) G A([, then x(0+) G S - S~ - A~; if x(0) G AZ, then x{0+) G S~ ; if 
x(0) G Aq , then x(0+) G S - §+ - A+. 
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Proof. We only treat the first two cases, because the reasoning for the 
last two is the same. If x(0) G A+, then d'(x(0)) > by Lemma 7.1, which 
implies the result. If x(0) G Aq , then d'(x(0)) = by Lemma 7.1. To see 
why we cannot have x(0+) G A - U note that then tti^x) would jump 
from 1 to 0, which would cause a jump in d'(x) because of Lemma 7.1 and 
the inequality in (6.13). □ 

We now are ready to establish local Lipschitz continuity. 

Definition 7.1. (local Lipschitz continuity) A function f : Q.2 — > K m , 
where Q± C O2 f= is locally Lipschitz continuous on Qi within Q2 if, f or 
every «q G Oi, there exists a neighborhood U C fi 2 °f v such that f restricted 
to U is Lipschitz continuous; i.e., there exists a constant K = K(U) such 
that — 7(^2)11 < K\\vi — V2W for every v±,V2 G U. 

Theorem 7.1. The function \& in (4.2) is locally Lipschitz continuous 
on § + within S + , on §~ within S~, on A within S b , on A + within S — S~ 
and on A~ within S — S + . 

Note that all of S is covered by the five cases in Theorem 7.1. Also note 
that, in each case, when we conclude that '3/ is Lipschitz continuous on fli 
within 0,2, if a solution exists starting at some point in S7i , then it will 
necessarily remain in ^2 for a short interval, by the reasoning above. We 
postpone the relatively long proof until §12. 

7.2. Global Existence and Uniqueness. This section is devoted to com- 
pleting the proof of Theorem 5.2 by establishing global existence and unique- 
ness. We first observe that, in general, one overall differentiable solution to 
the ODE over [0, 00) may not exist. From either S~ or § + , the solution x 
can hit S b , i.e., one of the sets A, A~ and A + , and have drifts that are in- 
consistent with the drifts in the new destination set. For example, in S + we 
necessarily have 7Ti 2(3?) = 1- However, in general there is nothing preventing 
x(t) — > xfo), where x(t) G S + with iri^x^t)) = 1 but also 5+(x(t)) < < 
S-(x(t)) while x(tb) G A, necessarily with 5+(x(tb) < < £_ (a: (£&)). The 
probability 7Ti ) 2(x(t)) jumps instantaneously from 1 to some value strictly 
between and 1 when A is hit. A numerical example is given in the ap- 
pendix. To treat that case, We can start a new ODE at this hitting time of 
A. 

We first show that the possible values of x are contained in a compact 
subset of S, provided that the initial values of the queue lengths are con- 
strained. That is accomplished by proving that a solution to the IVP (4.3) 
is bounded. We use the notation: a V b = max{a, 6}. 
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Theorem 7.2. (boundedness) Every solution to the IVP (4.3) is bounded. 
In particular, the following upper bounds for the fluid queues hold: 

(7.5) qi (t) < qf = ft(0) V Xi/Oi t>0, i = l,2. 

Proof. Since < z\$, < mi and qi > in S, we only need to establish 
(7.5). To do so, it suffices to consider the bounding function describing the 
queue-length process of each queue in a modified system with no service 
processes, so that all the fluid output is due to abandonment, which produces 
a simple one-dimensional ODE for each queue; for the remaining details, see 
§D in the appendix. □ 

PROOF of Theorem 5.2 (ii). It follows from Theorem 5.2 (i) established 
above, and Theorems 7.1 and 7.2, that any solution x on [0,5) can be 
extended to an interval [0,(5'), 8' > 8 (even 5' = oo), with the solution 
{x(t) : t G [0,5')} again being unique, provided that that the solution x 
makes no transitions from §> — S & to A, causing a discontinuity in ~ki^(x) 
and thus in (4.2). (See Theorem 3.3 in [9] and its proof for supporting 
details.) 

Moreover, the solution in S + or §>~ has a left limit at the time it hits 
A. The left limit exists because, by Theorem 7.2, the solution is bounded, 
and because the derivative in either S + or §~ is bounded, by (4.2). At each 
such hitting time, a new ODE is constructed starting in A. That ensures the 
overall continuity of x. In general, there can be accumulation points of such 
hitting times of the set A from 8 — S b . However, any such accumulation point 
t must be in either A + or A - . That is so, because there then are sequences 
{4 : n > 1}, i = 1, 2 with x(t l n ) G § - S b and x{t 2 n ) G A for all n with t l n f t 
and x(t l n ) — > x(t) G § b as n — > oo for i = 1,2. Finally, by Theorem 7.1, the 
function \& is locally Lipschitz continuous at each point in A + U A~. Hence, 
the solution x must actually be differentiable at each of these accumulation 
times of hitting times. As a consequence, x is continuous and differentiable 
almost everywhere throughout [0, oo). □ 

In the proof of Theorem 5.2 (ii), just completed, we have also established 
the following result. 

Theorem 7.3. (extension to a global solution) Let x be the unique dif- 
ferentiable solution to the IVP (4.3) on an interval [0,8), established in §7.1. 
If it is known that the solution can never transition from S + or E>~ to A then 
there exists a unique differentiable solution to the IVP (4.3) on [0, oo). 

8. Fluid Stationarity. We now define a stationary point for an ODE 
and then show that there exists a unique one for the ODE (4.2). We then 
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give conditions under which the fluid solution x = {x(t) : t > 0} converges 
to stationarity as t — > oo. In §10, we show that it does so exponentially fast. 

Definition 8.1. (stationary point for the fluid) We say that x* is a 
stationary point for the ODE (or fluid model) if x(t) = x* for allt > when 
x(0) = x* . That is, x* is a stationary point if^(x*) = for ^ in (4.1) and 
(4.2). If x(t) = x* for all t, then we say that the fluid solution is stationary, 
or in steady state. 

8.1. Characterization of the Stationary Point. By definition, a station- 
ary point x* G S satisfies ^/(x*) = 0. From (4.2), we see that this gives 
a system of three equations with three unknowns, namely, qf, q 2 and z* 2 . 
The apparent fourth variable 7rJ 2 = ^i,2( x *) is a function of the other three 
variables and its value is determined by x* . In principle, the three equations 
in ty(x) = can be solved directly to find all the roots of ^> . However, 7r£ 2 is 
a complicated function of x* having the complicated closed-form expression 
in (6.14) and (6.17). 

Theorem 8.1 below states that, if there exists a stationary point for the 
fluid ODE (4.2), then this point is unique, and must have the specified 
form. The uniqueness of x* is proved by treating fourth variable, 

and adding a fourth equation to the three equations *$>(x) = 0. However, it 
does not prove that a stationary point exists. In general, the solution 7r| 2 
we get from the system of four equations may not equal to TTi t2 (x*), for the 
function -k\^ defined in (3.7). The existence of a stationary point is proved 
in the next section. 

The proof of existence is immediate from the proof of uniqueness when 
7Ti 2(2;*) is known in advance to be or 1, with the value determined. That 
occurs everywhere except the region A; it occurs in the two regions S + and 
S~, but it also occurs in § b — A. Since the QBD is not positive recurrent in 
S b — A, it follows that vri 5 2(x*) can only assume one of the values, or 1, 
achieving the same value as in the neighboring region S + or §>~. (We omit 
detailed demonstration.) But we will have to work harder in A. 

We now focus on uniqueness. Although tv\ 2 is treated as a variable, we 
still impose conditions on it so that it can be a legitimate solution to (3.7). 
In particular, if q\ — rq\ > n then we let it\ 2 = M ^ Q* ~ r l2 < K ' then 
we let 7r* 2 = 0. Equation (8.3) below shows that < tt\ 2 < 1 whenever 
q* — rq 2 = k, i.e., whenever x* G S b . 

For a, b G R, recall that a V b = max{a, 6} and let a A b = min{a, b}. Let 

,g -g z _ 6>2(Ai - mifii y i) - r6i(X 2 - 7712^2,2) - O1O2K 

?~#l/i2,2 + #2/^1,2 
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Theorem 8.1. (uniqueness of the stationary point) There can be at most 
one stationary point x* = (q*, q%, z\ 2 ) for the IVP (4.3), which must take 
the form 

(8.2) 

z 12 = 0VzAm 2 , q 1 = , q 2 ~ 



for z in (8.1). Moreover, 
(8.3) ttU 



d 1 ' ^ 2 



^1,2^1,2 



Ml,2^*,2 + ( m 2 - 4,2)/^2,2 



Proof. We start with (8.3). This expression is easily derived from the 
third equation in (4.2), by equating z\^ 2 {t) to zero. Observe that if z* 2 = m 2 
then tt^ 2 in (8.3) is equal to 1, and if z* 2 = then 7r* 2 = 0. Now, by 
plugging the value of it\ 2 in the ODE's for <?i(i) and c/2(£) in (4.2) we get 
the expressions of q* and q\ in (8.2). We now have the two equations for 
the stationary queues, but there are three unknowns: z* 2 , q\ and q 2 . We 
introduce a third equation to resolve this difficulty 

Consider the following three equations with the three unknowns: z,q±(z) 
and q 2 {z). (here q\ and q 2 are treated as functions of the variable z, not to 
be confused with the fluid solution which is a function of the time argument 
t.) 

,s Ai - (tii,imi - m_ 2Z f >. A 2 - M2,2(m 2 - z) 

(8.4) giW = 0, ' q2{z) = J 2 ' 

K = qi{z) - rq 2 (z). 

Notice that qi(z) is decreasing with z, whereas q 2 (z) is increasing with z. 
Thus, there exists a unique solution to these three equations, which has z 
as in (8.1). We can recover x* from the solution to (8.4), and by doing so 
show that x* is unique and is always in one of the three regions §~, S + or 
S b (so that x* £ §>). 

Let (qi(z), q 2 (z), z) be the unique solution to (8.4). First assume that 
z > 777-2, which implies that q 2 (z) > 0, and, by the third equation, qi(z) > K. 
By replacing z with m 2 , q±(-) is increased and q 2 (-) is decreased (but is 
still positive), so that qi(m 2 ) — rq 2 (m 2 ) > k (and, trivially, qi(m 2 ) > k, 
92(^2) > 0). This implies that x* = (qi(m 2 ), q 2 (m 2 ),m 2 ) £ E + and, if it is 
indeed a solution to *$>(x) = 0, then x* is the unique stationary point for the 
ODE. 
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Now assume that the unique solution to (8.4) has z < 0. By replacing z 
with we have qi(0) < qi(z) and ^(O) > q2 (z), which imply that qi(0) — 
r ?2(0) < k. Now, since qi(0) = q^ we have that qi(0) > k by Assumption A. 
This implies that qi(z) > k, which further implies that rq2(z) = qi(z) — k > 
0, so that rq 2 (0) > rq 2 (z) > 0. Taking x* = (<?i(0), g 2 (0), 0), we see that 
x* £ and if x* is indeed a solution to *$>(x) = 0, then x* is the unique 
stationary point for the ODE. 

Finally, assume that the solution x(z) = (qi(z),q 2 (z), z) to (8.4) has < 
z < m 2 . To conclude that x(z) is in S 6 we need to show that q(z),q 2 (z) > 0, 
so that q\ = qi(z) and q 2 = q 2 {z) are legitimate queue-length solutions. We 
now show that is the case under Assumption A. 

Let S^ = m 2 - A 2 //i 2 ,2- Note that, if 5f > 0, then 5f = a£, for in (5.1). 
We start by rewriting qi(z) and q 2 (z) in (8.4) as 

(8.5) q 1 (z)=q" 1 -^z, q 2 {z) = !f*(z - S a 2 ). 
Now, it follows from Assumption A that 

(8.6) K < ql - ^fs a 2 <qt- ^S a 2 , 

where the second inequality follows trivially, since S 2 < s 2 . From the third 
equation of (8.4), k = qi(z) — rq 2 (z). Combining this with (8.5), we see that 

(8.7) k = qi (z) - rq 2 {z) = qf - ^z - r^(z - S a 2 ). 
Combining (8.6) and (8.7), we get 

Ql - r ^^\ z - b 2) ^ 1l ^^2) 

Cl tl 2 (71 

which is equivalent to 

°-(^ + r lf) {Z - S ^ 

This, together with the fact that the solution has z > 0, implies that z > 
max{0, S 2 } = s 2 . It follows from (8.5) that q 2 (z) > and, by using the third 
equation in (8.4) again, qi(z) = rq 2 {z) + k > k > 0. □ 

An immediate consequence of the proof of Theorem 8.1 is that, in order to 
find the candidate stationary point x* , one has to solve the three equations 
in (8.4). The next corollary summarizes the values x* may take, depending 
on its region; the proof appears in the appendix. 
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Corollary 8.1. Let x* = (qi,q2,Zi 2 ) be the point defined in Theorem 
8.1. 

1. If x* G S b , then, for z defined in (8.1), 

f 1^2(9? - «) - ^6 | i(A 2 - ^2,2^2) 



Zl ' 2 ^ r0l/X 2 ,2 + #2^1,2 

6>ie 2 (g?-rgf-K) 
r6llJ,2,2+02Hl,2 ' 



*/ <z 2 a > 0, 4 = 0. 



6>ie2K+r'M2.2sg/e 2 -«) ., a _ n „a ^ n 
r0i M2 ,2+02Ml,2 ' 7 92 ~~ ' 2 > 

Ai - mi/in - 4 2 ^ii 2 A 2 - (m 2 - z\ 2 )/x 2 ,2 

9l = a ■ > 12 = n 1 • 



2. Ifx* = §+, then 



* * Ai - mi/ii,i - m2/ii,2 * A 2 

«1,2 = m 2, gi = 7 , 9 2 = 

#1 #2 

5. J/x* e §~, tfien 

n * Ai - min^i „ A 2 - m 2 //2,2 
^1,2 = 0, ?i = , q 2 = 2 • 

If x* £ S + , as in (ii), then the system does not have enough service 
capacity to keep the weighted difference between the two queues at k, even 
when all agents are working with class 1. In this case, the only output from 
queue 2 is due to abandonment, since no class-2 fluid is being served (in 
steady state). Queue 2 is then equivalent to the fluid approximation for an 
M/M/oo system with service rate 82 and arrival rate A2. On the other hand, 
queue 1 is equivalent to an overloaded inverted-!^ model: a system in which 
one class, having one queue, is served by two different service pools. 

The next corollary gives necessary and sufficient conditions for x* to be 
in each region. It shows that the region of x* can be determined from rate 
considerations alone. We give the proof in the appendix. 

Corollary 8.2. Let x* be as in (8.2). Then 
1. x* £ S 6 if and only if 



12 vi 



x* 6 A if and only if both inequalities are strict. 



20 



O. PERRY AND W. WHITT 



2. x* G §+ if and only if qf-n> r -^ + ^f^. 

3. x* G S _ if and only if rq^ > qf — k. 

Remark 8.1. (most likely region in applications) It follows from Corol- 
lary 8.2 that, in applications, A is the most likely region for the stationary 
point when the system is overloaded, provided that the arrival rates are 
about 10 — 50% larger than planned during an overload incident. Typi- 
cally, a much higher overload is needed in order for the stationary point 
to be in S + . As an example, consider the canonical example from [15]: 
There are 100 servers in each pool, serving their own class at rates fin = 
^2,2 = 1- Type-2 servers serve class- 1 customers at rate /ii^ = 0.8. Also, 
0\ = 02 = 0.3, r = 0.8 and k = 0. Suppose that class 2 is not overloaded 
with A2 = 90. Then, for the stationary point to be in S + , we need to have 
Ai > f*i imi + Mi,2^2 + 6^X2/62 = 252, i.e., the class-1 arrival rate is 252% 
larger than the total service rate of pool 1. If A2 > 90, especially if pool 2 is 
also overloaded, then Ai needs to be even larger than that. 

8.2. Existence of a Stationary Point. We have just established unique- 
ness of the stationary point in S, and characterized it. In the process, we 
have also established existence in § — A. Now we will establish existence of 
the stationary point in A. First, we calculate the drift rates at x* G A. 

Lemma 8.1. (the drift rates at x*) For x* in Corollary 8.1 (i), where 

< Z* >2 < 1712, 

(8.9) 5+(x*) = -(i+A0/i2,2(m 2 -z 1)2 ) < °> d -( x *) = +(j+^Vi,2< 2 > 0. 
PROOF. Substitute x* in Corollary 8.1 (i) into (6.9), using (6.1)-(6.4). □ 
We now are ready to prove existence. 

Theorem 8.2. (existence) If the model parameters produce x* G A, i.e., 
as in Corollary 8.1 (i), where < z\ 2 < wi2, then x* is the unique stationary 
point. 

Proof. We will prove that there must exist at least one stationary point. 
Given that result, by Theorem 8.1 and Corollary 8.1, there must be exactly 
one fixed point and that must be the x* given there. To establish existence, 
we will apply the Brouwer fixed point theorem. It concludes that a continu- 
ous function mapping a convex compact subset of Euclidean space ffi fc into 
itself has at least one fixed point. We will let our domain be the set 

(8.10) C(r]) = {x G AnB : 5+ (x) < -?? and 6-(x) > i]} 
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for an appropriate small positive 77, where B = [0, q\ d ] x [0, q 2 d ] x [0, m 2 ] with 
qbd being the bound on % from Theorem 7.2. Choose t] sufficiently small 
that x* G C(r]); that is easily ensured by Lemma 8.1. Since the rates in 
(6.1)-(6.4) and the drift in (6.9) are linear functions of x, we see that C(ij) 
is a convex subset of A for each r\ > 0. Since the inequalities in (8.10) are 
weak, the set is closed. The intersection with B guarantees that the set C(r]) 
is also bounded. Hence, C(jj) is compact. 

By Theorem 5.2, for any x(0) G C(rf), there exists a unique solution to the 
ODE over [0, 5] for some positive 5. Hence, for any t with < t < 5, the map 
from x(0) to x(t) is continuous; see §2.4 of [19]. Let x* L = (q^ L , q 2 L , z* 2 L ) 
and x* v = (q{u,q* 2iU ,zl 2>u ), where q* 1L = qf-e, q* 2 L = q* 2 -e, z* 2L = zfj-e, 

1*,u = Qi + ^ lt,u = 1*2 + e and z t,2,u = 4,2 + e - Let <k c (v) ->■ C{rf) be 
the continuous function defined by 4>t(x(0)) = (qi,t, Q2,t, z i,2,t)i where 

(8.11) qi^t = Qi(t) V q*^ A q*y and z\,2,t = z l,2,t V Zi i2 ,l ^ z l,2,Ui 

for i = 1,2. We can choose i] > and e > sufficiently small so that, first, 
x* G C{rf) and, second, that Xi t t G C{rj) for each x(0) £ C(rj). Hence, the 
pair {C{rf),4>t) satisfies the conditions for the Brouwer fixed point theorem. 
Hence, there exists x(0) G C(rj) such that x(t) = x(0). 

Now let {t n : n > 1} be a sequence of time points decreasing toward 0. 
We can apply the argument above to deduce that, for each n, there exists 
x n (0) in C(rj) such that x n (t n ) = x n (0). However, from the ODE, we have 
the relation \x(t) — x(0) — ^(x(0))t\ < Mt 2 for all sufficiently small t. Since 
{x n (0) : n > 1} is bounded, there exists a convergent subsequence. Let x(0) 
be the limit of that convergent subsequence. For that limit, we necessarily 
have ^(x(0)) = 0. Hence, that x(0) must be a stationary point for the ODE. 
By Theorem 8.1, we must have x(0) = x* . □ 

8.3. Global Asymptotic Stability. Having a unique stationary point does 
not imply that a fluid solution necessarily converges to that point as t — > 00. 
It does not even guarantee that a solution to the IVP (4.3) is asymptotically 
stable in the sense that, if ||x(0) — x*\\ < e, then x(t) — > x* as t — > 00, no 
matter how small e is. In fact, there is not even a guarantee that x(t) will 
remain in the e-neighborhood of x* for all t > 0. We will establish all of 
these properties in Theorem 8.3 below by showing that x* in §8.1 is globally 
asymptotically stable, as defined below: 

Definition 8.2. (global asymptotic stability) A point x* is said to be 
globally asymptotically stable if it is a stationary point and if, for any initial 
state x(0) and any e > 0, there exists a time T = T(x(0),e) > such that 
\\x(t) -x*\\ < e for allt> T. 
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Global asymptotic stability goes beyond simple convergence by also re- 
quiring that the limit be a stationary point. 

Theorem 8.3. (global asymptotic stability of x*) The unique stationary 
point x* is globally asymptotically stable. 

The proof of Theorem 8.3 relies on Lyapunov stability theory for deter- 
ministic dynamical systems; see Chapter 4 of Khalil [9]. Let E be an open 
and connected subset of W 1 containing the origin. We use standard vector no- 
tation to denote the inner product of vectors a,b € M. n , i.e., a-b = Y27=i a i^i- 

Definition 8.3. (Lie derivative) For a continuously differentiable func- 
tion V : E — > R, and a function ^ : E — > W 1 , the Lie derivative of V along 
^ is defined by 

dV 

V(x) = -z-V(x) = VV • 
ox 

where W = (Jj-, • • • , is the gradient of V . 

Definition 8.4. (Lyapunov-function candidate) A continuously differ- 
entiable function V : E — )■ R is a Lyapunov-function candidate if: 

1. V(0) = 

2. V{x) > for all x in E - {0} 

In proving Theorem 8.3 we use the following theorem, which is Theorem 
4.2 pg. 124 in [9]: 

Theorem 8.4. (global asymptotic stability for nonlinear ODE) Let x = 
be a stationary point of x = ^(x), ^ : E — > R n , and let V : R" — > R be a 
Lyapunov-function candidate. If 

1. V(x) — > oo as \\x\\ — > oo and 

2. V(x) < for all x^O, 

then x = is globally asymptotically stable as in Definition 8.2. 

Notice that, under the conditions of Theorem 8.4, the Lyapunov-function 
candidate V provides a form of monotonicity: We necessarily have V(0) = 
and V(x(t)) strictly decreasing in t for x(t) ^ 0. To elaborate, we introduce 
the notion of a y-ball. We say that /3y(a) is the a V-ball with center at x* 
and radius a if 

(8.12) p v (a) = {x£R n : \\V(x) - V(x*)\\ < a}. 

If x(to) S /3y(a) f° r some a > and to > 0, then x(t) 6 /3y(a) for all t > to- 
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Proof of Theorem 8.3. Theorem 8.4 applies directly only within one 
region, starting at a point in § + , S~, A, A~ or A + . However, we will show 
that the same Lyapunov function V can be used in all regions, leading to 
global decrease of V as x* is being approached. 

Let x be the unique solution to (4.3). Let x* = (</*, 92, 2*2) be the sta- 
tionary point for the system (4.1). Without loss of generality, we perform a 
change of variables and define a new system whose unique stationary point 
is x = 0. To this end, let y = x — x* so that y = x = ^Sf(x). Hence, 
q>(x) = fy(y + x*) = g(y) and we have that g(0) = ^(0 + x*) = ^(x*) = 0. 
That is, if x* is a stationary point for the original system x = ^(x), then 
the stationary point for the new system, y = g(y), is y* = 0. We distinguish 
between two cases: (i) W,2 > ^2,2 and {%%) pL\ )2 < \i 2 ^ 2 . 

(i) First, if \i\ >2 > ^2,2, then choose V\{x) = x\ + X2 and apply its Lie 
derivative along g (y) = ^f(y+x*) where y+x* = (gi (0+9*> 92(0+92' z i,2(0 + 
z\ 2 ) an d x* is given in (8.2). By the definition of the Lie derivative, V\{y) 
is equal to the inner product 

Vi(y) = (1,1,0) • (gi(0,9 2 (0,M0)' =9i(0 + 9 2 (0, 

for qi, q 2 and Zi j2 in (4.2), after the change of variables. Let £1,2 (0 = 21,2(0 + 
z*. Then, for x* = (9^,92,^1,2) as in ( 8 -2) 

V\{y) = Ai - mi^i - vri i2 (y(t))[zi >2 (0^i,2 + (m 2 - 21,2(0)^2,2] 

- 0i(9i(O + ?*) - (1 - vri, 2 (y(t)))[(m 2 - 2i, 2 (0)M2,2 + ^(O/M 
+ A 2 -02(92(0 + 9*) 

= Ai + A 2 - mi/^1,1 - m 2/ u 2 ,2 + 21,2(0^2,2 + 2*^2,2 - 21,2(^)^1,2 

- 2^1,2 - #i9i(0 - #i9i - #292(0 - #292 
= -6>i<?i(0 - 6»29 2 (0 - ^1,2(0(^1,2 - ^2,2)- 

Thus, Vi(y) < for all t/£t 3 unless y = 0. 

(ii) When /ii,2 < ^2,2, there exists a B > 1 such that /i 2t2 = B/ii^. 
We next show that for any C > B the candidate-function V 2 (x) = Cx\ + 
x 2 + (C — 1)^3 is a Lyapunov function. The Lie derivative of V 2 {x) for the 
modified system g(y) is 

V 2 (y) = (C, 1, C - 1) • (91 (0, 92(0, 2i )2 (0) = Cgi(t) + g 2 (t) + (C - l)ii, 2 (0- 
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Hence, 

V 2 (y) = C[Xi- - vri i2 (y(t))(zi i2 (t)Mi,2 + {m 2 - 5i j2 (i))// 2j2 )] 

- &l{<ll(t) + q\) + A 2 - e 2 {q 2 {t) + ^) 

- (1 - 7Tl,2(y(i)))(5l,2(*)Ml,2 + (m 2 - Sl, 2 (t)M2,2)) 

+ (C - 1) [7ri )2 (i/(t))(m 2 - 2i,2(0)/i2,2 - (1 - 7ri, 2 (y(t)))zi, 2 (i)/ii )2 ] 
= -C0igi(t) - 2 g 2 (t) - z li2 (t)(Cfn, 2 - m 2j2 ), 

so that V^y) < for all y / 0. 

By Theorem 8.4, y* = is globally asymptotically stable for the modified 
system g{y). Hence, x* is globally asymptotically stable for the original 
system That is, for every initial value x(0) we have that x(t) — > x* . □ 

8.4. Staying in 8. We also use the Lyapunov argument to prove Theorem 
5.1, i.e., show that the solution to the ODE can never leave S. 

PROOF of Theorem 5.1. We are given x(0) G S. Consider t > 0. It is easy 
to see that, if z\ t2 {t) = 0, then ii, 2 (t) > 0, so that z\ t2 (t+) > 0. Similarly, if 
z i,2(t) = m 2 , then ii )2 (t) < 0, so that ^i i2 (t+) < m 2 . 

Turning to the queues, note that to leave § at time t+ we must have 
qi(t) = k or q 2 (t) = (or both). If qi(t) = k and q 2 {t) > 0, then x(t) G S~ 
so that TTi, 2 (x(t)) = 0. Plugging this value of 7Ti i2 (x(t)) in the ODE for q\{t) 
in (4.2), we see that q\{t) > Ai — ii\\m\ — 0\k > by Assumption A. 
Hence, qi(t) is nondecr easing. If q±(t) > k and q 2 (t) = 0, then x(t) £ S + and 
7r i,2( 2; (i)) = 1) which gives q 2 (t) = A 2 > 0. Hence q 2 is increasing at time t. 

Now consider the case qi(t) = k and q 2 (t) = 0, so that x(t) G E> b . For one 
of the queues to become negative at time t+, we need to have its derivative 
be negative at time t. We will consider various subcases. 

First assume that q\(t) < and q 2 (t) > 0. In that case qi(t+) < q 2 (t+), 
so that 7Ti, 2 (x(t+)) = 0. Plugging this value of 7ri ]2 (x(i+)) in the ODE 
(4.2), together with qi(t+) = k, we see that qi(t+) > by Assumption A. 
Next assume that q\(t) > and q 2 (t) < 0. Then qi(t+) > q 2 (t+), so that 
7Ti )2 (x(t+)) = 1. Plugging this value of vri j2 (x(t+)), together with q 2 {t+) = 
0, we see that q 2 {t+) > 0. 

We finally consider the remaining more challenging subcase: qi(t) < 
and q 2 (t) < 0. We will show that this subcase is not possible. To that end, 
we further divide this case into three subcases: x(t) G A + , x(t) G A~ and 
x(t) G A. (Recall that = A U A + U A~.) However, x(t) cannot be in A~, 
since then 7Ti >2 (x(i)) = 0, which implies that qi(t) is nondecreasing (plug 
^i,2(x(t)) = and qi(t) = k into the ODE (4.2)). Moreover, x(t) cannot 
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be in A + , since then ni t 2(x(t)) = 1, which implies that q2(t) is strictly 
increasing. 

Now assume the remaining possibility, x(t) E A, and recall that is 
Lipschitz continuous in A, so that the Lyapunov argument holds over [t, t + 
rj), for some rj > 0. Specifically, the Lyapunov function V is monotone 
increasing in x(t), because x* > 0. (The inequality holds componentwise.) 
If Ml,2 > M2,2> then we take the Lyapunov function Vi(x(t)) = qi(t) + q 2 (t). 
The monotonicity of V\ at x(t) implies that at least one of the queues must 
be increasing, which contradicts the assumption that the derivative of both 
queues is negative at t. If //x,2 < ^2,2; then we take the Lyapunov function 
V 2 (x(t)) = Cqi{t) + q 2 (t) + {C - ]>i, 2 (X). We then choose C = 1 + e with 
e small enough, such that V2(x(t)) < (assuming the derivatives of both 
queues are strictly negative at t). Once again, this contradicts the positive 
monotonicity of V at x(t). This concludes the proof. □ 

9. Exponential Stability. 

Definition 9.1. (exponential stability) A stationary point x* is said 
to be exponentially stable if there exist two real constants (3 > such that 

\\x{t) - x*\\ < 0||a?(O) - x*\\e-P\ 

for all t > and for all x(0), where \\ ■ \\ is a norm on W 1 . 

We use Theorem 3.4 on p. 82 of Marquez [11], stated below. 

Theorem 9.1. (exponential stability of the origin) Suppose that all the 
conditions of Theorem 8.4 are satisfied. In addition, assume that there exist 
positive constants K\, K 2 , K3 and p such that 

Ki\\x\\ p < V(x), <K 2 \\x\\ p and V(x) < -K 3 \\x\\ p . 

Then the origin is exponentially stable, and 

\\x(t)\\ < \\x(0)\\ {Ki/Kxf'Pe-^/ 2 ^ 1 for all t and x(0). 

We use the L\ norm: ||x|| = \x\\ + \x 2 \ + \x%\ for x G M 3 . 

Theorem 9.2. (exponential stability of x*) Each x* in S is exponentially 
stable. 

1. If Hx,2 > fJ-2,2, then 

\\x(t)-x*\\<\\x(0)-x*\\e- ( - K3/2)t for all t and x(0) 

for all x(0) £ S and t > 0, where K3 = max{^i, 9 2 , /ii,2 — ^2,2}- 
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2- If H2,2 = Bn^ 2 , B>1, then 

\\x{t) -x*\\< ||x(0) - xKC/A'Oe-^/ 2 )' 

for all x(0) £ S, t > and C > B, where K\ = min{l,C — 1} and 
K A = m&x{C9i, 6 2 , (C/ii )2 - ^2,2)} ■ 

Proof. As in the proof of Theorem 8.3, Theorem 9.2 applies directly 
only within one region, starting at a point in S + , §>~, A, A~ or A + . However, 
again, the same Lyapunov function V can be used in all regions. 

We consider the two cases in turn: (i) In the proof of Theorem 8.4, V\(x) = 
x\ + X2 was shown to be a Lyapunov function with a strictly negative Lie 
derivative. Since x > 0, we can take K\ = K2 = 1 and p = 1. Since 
V\(x) = -8iqi(t) - 9 2 q2{t) - {^1,2 - ^2,2)^1, 2{t), we can take K 3 in (??), and 
the result follows from Theorem 9.1. 

(ii) We use the Lyapunov function V2(x) = Cx\ + X2 + (C — l)x 3 . Then 
Ai||x|| < V2(x) < C||x|| for K\ = min{l,C — 1}. From the proof of Theorem 
8.4, we know that V 2 (x) = -C6iqi(t) - 02?2(i) - (C/"i,2 - ^2,2)^1,2(^)1 so 
that V 2 (x) < -A 4 ||x||. □ 

10. Conditions for State-Space Collapse. In this section we give 
ways of verifying that x lies entirely in A, given that x(0) and x* are both 
in A. In the appendix we provide conditions for the solution to eventually 
reach A after an initial transient. The results here are intended to apply 
after this initial transient period has concluded. 

Theorem 10.1. (sufficient conditions for global SSC) Let v = ^^^^2,2, 
and suppose that x(0) £ A. Also assume that 

(10.1) q 2 (0) < X2/O2 and 9l (0) < (Ai - m^i,i)/fli. 

If, in addition, the following inequalities are satisfied, then the solution to 
the IVP (4.3) is in A for all t: 

(10.2) (i) \\ < vrri2 + mifii t i and (ii) A2 > vm,2- 

Proof. We start by showing, under Condition (i), that 5+(x(t)) in (6.9) 
is strictly negative for each t. For a fixed t, 

S + (x(t)) = j (A?(t) - ${t)) + k (A«(t) - M «(t)) < 

if and only if 
(10.3) 

(M2,2 - ^1,2)^1,2(0 - "12^2,2 < -(Ai - mi/ii,i) + r(A 2 - 9 2 q2{t)) + 6iqi(t). 
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If A*2,2 > Mi,2j then the left-hand side (LHS) of (10.3) is maximized at 
z i,2(t) = m2, and is equal to — [ix^m?.. If A*2,2 < Mi,2j the the LHS is maxi- 
mized at 2i,2(i) = 0, and is equal to —fi 2:2 m 2 . When [i 22 = Mi,2 the LHS is 
equal to — /« 2i2 m 2 = — ^1,2^2 • Overall, the LHS of (10.3) is smaller than or 
equal to —vm 2 . 

Since 92(0) < A 2 /# 2 , we conclude, using the bound in (7.5), that 9 2 q 2 {t) < 
A2 for all t > 0. This, together with the fact that qi(t) > for all t, implies 
that the RHS of (10.3) is larger than or equal to — (Ai — mi/i^i), so that 

(a*2,2 - ^1,2)^1,2^) - M2,2"T-2 < -vm 2 < -(Ai - m 1 fi 1) i) 

< -(Ai - mi/ii,i) + r(A 2 - ftjga(t)) + 0igi(t) 

where the second inequality is due to condition (i). 

To show that condition (ii) is sufficient to have 6-(x(t)) > for all t, fix 
t > and note that, for <5_(x(i)) in (6.9), we have 

S-(x(t)) = j (A^(t) - M ^(t)) + fc (A W (t) - > 

if and only if 
(10.4) 

r(m,2 - ^2,2)^1,2(0 + rn2,2 m 2 > -(Ai - mi/ii,i) + r(A 2 - 2 q 2 (t)) + 0\qi(t). 

It is easy to see that the LHS of (10.4) has a minimum value of r(fj,± 2 A 
^2,2)^2 = rum 2 . By essentially the same arguments as in Theorem 7.2 
we can show that <Zi(i) < <7i(0) V (Ai — m\^i t i) /0\. Since we assume that 
<7i(0) < (Ai — mifii t i)/8i, we have the bound gi(t) < (Ai — mifj,i t i)/9i for 
all t > 0. With this bound, we see that the RHS of (10.4) is smaller than or 
equal to rA 2 . Overall, we have 

r(m,2 ~ ^2,2)^1,2(0 + rH2,2m 2 > rvm 2 > r\ 2 

> -(Ai - mi/ii,i) + r ( A 2 - + 0iqi{t), 

where the second inequality is due to Condition (ii). Since (6.10) holds for 
all t > 0, we also have < iri )2 (t) < 1 for all t. Hence, every solution to the 
IVP in (4.3) must lie entirely in A. □ 

For x* £ A, we will now show that there exist a > and T = T(q), such 
that global SSC can be inferred once ||x(T) — sc*|| < a. We exploit the drift 
rates at stationarity, defined by 5* + = S+(x*) and 5*L = <5_(a;*). It follows 
from the expressions in (6.9) that 
(10.5) 

6* + = 5 + {x*) = -/i 2 , 2 (r + l)(m 2 - z* h2 ), 6*_ = 5_(x*) = Mlj2 (r + 1)< 2 . 
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Thus, if < z\ 2 < m 2, then the positive recurrence condition (6.10) holds 
at the stationary point x* . (This agrees with (8.3) which has < n\ 2 < 1 if 
and only if < z* 2 < rri2-) 

In the next theorem we give explicit expressions for a. For reasonable 
rates, such as will hold in applications, a is quite large. In fact, in the nu- 
merical example considered in §11.3 we show that, typically in applications, 
a is so large, that we can infer that x lies entirely in A without even solving 
the IVP; i.e., x(0) G Pv{a). 

Theorem 10.2. Suppose that x* G A and let £ = min{|#!j_|, SI.}. 

1- If M22 > Ml 2> then let a = £/V#2 

2- If 1^2,2 < Ml, 2, then let a = u>/iere ? = ^ 1)2 - M2,2 + #1 + r#2 > 0. 

In both cases, if there exists T > suc/i that x(T) G /3y(a), i/ien {x(t) : t > 
T} lies entirely in A. 

Proof. We use /3y(a), the a V-b&ll with center at x* and radius a, 
in (8.12). To find a proper a for the I^-ball /3y(a), we once again use the 
conditions (10.3) and (10.4). We first show how to find a for the case //2,2 = 
B/ii t 2 for some B > 1, i.e., when /ii^ < ^2,2- Recall (proof of Theorem 8.3) 
that in this case, V 2 (x) = Cxi + X2 + (C — l)x3 is a Lyapunov function 
for any C > B. Also, the Lyapunov function was defined for the modified 
system in which the origin was the stationary point. 

Let x* = {q\-,q 2l z\ 2 ) be the stationary point in A. First assume that, at 
some time T, V 2 (x{T)) = e x , i.e., Cq x {T) + q 2 (T) + (C - l)z ly2 (T) = e x . If 
x(t) G /3y 2 (ei) for all t > T, then it must hold that 

(10.6) 

Ql ~ -j^j < Ql(t) < Ql + -j^, 12 ~ £ i < Q2(t) < q 2 + ei and 

z *< 2 ~ cT~[ < Zl > 2 ® < z *^ 2 + cTTT' t - T - 

To make sure <5+(x(i)) < 0, we use (10.3), reorganizing the terms. We need 
to have 

(M2,2 - Mi,2)2i,2(*) + r9 2 q2(t) - 6\q x {t) < -(Ai - Mi.i^i) + r A 2 + H 2 , 2 m 2 . 
By (10.6), the above inequality holds if 

0*2,2 - Ml,2) (< 2 + + r0 2 (92 + e i) " #1 (9* " 

< - (Ai - /Ui^mi) + rA 2 + M2,2m 2 . 
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Plugging in the expressions for q\ , q\ and z* 2 , we see that we need to find 
an ei > such that 

(a*2,2 - 1*1,2) 1 + r0 2 ei + 61^ < H2,2{r + l)(m 2 - 2r* j2 ). 

We can take C as large as needed, so that the only term that matters on 
the LHS is r# 2 ei. Hence, we need to have 

At2,2(r + l)(m2-2 12 ) 1511 
ei < 



rt>2 

Similarly, to make sure that S-(x(t)) > 0, we use (10.4), reorganizing the 
terms. We need to have 

r(m,2 ~ 1*2,2) Zl,2(t) + r9 2 q2(t) - 6iqi{t) 
> -(Ai - m,imi) + r(A 2 - ^ 2 , 2 m 2 ). 

Using (10.6) again (with a different e 2 ), we see that it suffices to show that 

r(Mi, 2 - ^2,2) (z{ >2 + -^—^J + r9 2 {ql - e 2 ) - B\ (q\ + ^) 

> -(Ai - ii\,\m\) + r(A 2 - n 2 , 2 m 2 ). 

Once again, plugging in the values of q\ and z\ 2 , and taking C as large 
as needed, we can choose e 2 > such that 

r6> 2 r0 2 ' 

Hence, we can take a as in (i). 

For the second case, when [i\ %2 > /i 2j2 , we use the Lyapunov function 
V\{x) = x\ + x 2 . Using similar reasoning as above, we get 

H2, 2 (r + l)(m 2 - zU) \S*\ n 1>2 {r + l)z* 12 5*_ 

€\ < — = and e 2 < 



Mi,2 - M2,2 + 01 + ?~#2 C Ml,2 - A*2,2 + 6*1 + r6» 2 ? 

Hence, in this case we can take a in (ii). □ 

11. A Numerical Algorithm to Solve the I VP. 
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11.1. Computing tt\ i(x) at a point x. The QBD structure in §6.2 allows 
us to use established efficient numerical algorithms from [10] to solve for the 
steady state of the QBD to compute 771,2(2;), for any given x = x(t) £ A. 
We start by computing the rate matrix R = R(x). (To simplify notation, we 
drop the argument x, with the understanding that all matrices, are functions 
of x.) By Proposition 6.4.2 of [10], R is related to matrices G and U via 

(11.1) G = (-U)~ 1 A 2 , U = A 1 + A G and R = A (-U)~ 1 . 

In addition, the matrices G and R are the minimal nonnegative solutions to 
the quadratic matrix equations 

(11.2) A 2 + AiG + A G 2 = and A + RAi + R 2 A 2 = 0. 

Hence, if can compute the matrix G, then the rate matrix it! can be found 
via (11.1). Once R is known, we use (6.16) to compute ao- With ao and R 
in hand, 771,2(2;) is easily computed via (6.17). 

It remains to compute the matrix G. We use the logarithmic reduction 
(LR) algorithm in §8.4 of [10], modified to the continuous case, as in §8.7 
of [10]. The LR algorithm is quadratically convergent and is numerically 
well behaved. These two properties are important, because the matrix R(x) 
needs to be computed for many values of x when we numerically solve the 
IVP (4.3). From our experience with this algorithm, it takes fewer than ten 
iterations to achieve a 10 -6 precision (when calculating G). 

11.2. Computing the Solution x. To compute the solution x, we combine 
the forward Euler method for solving an ODE with the algorithm to solve 
for TTi t 2(x(t)) described above. Specifically we start with a specified initial 
value x(0), a step-size h and number of iterations n, such that nh = T. 
First, assume that 2i,i(0) = mi and 21,2(0) + £2,2(0) = m 2 , so that x(0) £ S. 
If 5(0) = (gi(0) - k) - rq 2 (0) > then n lj2 (x(0)) = 1. If 5(0) < then 
tti, 2(^(0)) = and if D(0) = then we check to see whether (6.10) holds. 
If it does, then x(0) £ A and we calculate 771,2 (x(0)) as described above. 
If x(0) £ S b - A then we can still determine the value of 771,2(^(0)) in the 
following way: If 5-(x(t)) = > 5+(x(t)), then we let 771,2 (a:(t)) = 0; if 
instead 5-(x(t)) > = 8+(x(t)), then we let 7Ti,2(x(i)) = 1. 

Given x(0) and 771,2(^(0)) we can calculate ^(x(0)) explicitly and perform 
the Euler step x(h) = x(0) + h^(x(0)). We then repeat the procedure for 
each k, < k < n — 1, i.e., 

(11.3) x((k + l)h) = x(kh) + h^(x(kh)), < k < n, 
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where x(kh) is given from the previous iteration, and ^f(x(kh)) can be com- 
puted once iri2(x(kh)) is found. 

If -2i,i(0) < m\ or ^1,2(0) + ^2,2(0) < ?7i2, so that x(0) ^ 8, we use the 
appropriate fluid model for the alternative region, as specified in the ap- 
pendix, where at each Euler step we check to see which fluid model should 
be applied. 

The forward Euler algorithm is known to have an error proportional to 
the step size h, and to be relatively numerically unstable at times, but it 
was found to be adequate. It would be easy to apply more sophisticated 
algorithms, such as general linear methods, which have a smaller error, and 
can be more numerically stable. The only adjustment required is to replace 
the Euler step in (11.3) by the alternative method. 

In the numerical example in §11.3 below we let the ratio be r = 0.8 = 4/5, 
so that all the matrices, used in the computations for n± 2, are of size 10 x 10. 
It took less than 10 seconds for the algorithm to terminate (using a relatively 
slow, 1 GB memory, laptop). The same example, but with r = 20/25, so 
that the matrices are now 50 x 50, the algorithm took less than a minute to 
terminate. Moreover, the answers to both trials were exactly the same, up 
to the 7th digit. In both cases, we performed 5000 Euler steps (each of size 
h = 0.01, so that the termination time is T = 50). It is easily seen that tti 2 
had to be calculated for over 4500 different points, starting at the time nip 
becomes positive (see Figure 2 in the following example). 

The validity of the solution can be verified by comparing it to simulation 
results, as in the example below and others in [15, 16]. There are two other 
ways to verify the validity: First, we can check that the solution converges 
to the stationary point x*, which can be computed explicitly using (8.2). 
Second, within A we can see that the two queues keep at the target ratio r, 
even though this relation between the two queues is not forced explicitly by 
the algorithm. 

11.3. A Numerical Example. We now provide a numerical example of the 
algorithm for solving the ODE in (4.1). In addition, we added the sample 
paths of the stochastic processes Q" an< ^ after scaling as in (3.2), on 
top of the trajectories of the solution to their fluid counterparts q\ and z\^. 

The model has the same target ratio r = 0.8 as in the example in §6.2 
with component rate matrices in (6.12). We chose a large queueing system 
with scaling factor n = 1000, so that the stochastic fluctuations do not to 
hide the general structure of the simulated sample paths. We let the ODE 
model parameters be mi = rri2 = 1, Ai = 1.3, A2 = 0.9, fi± : i = ^2,2 = 1, 
Mi,2 = M2,i = 0.8, 61 = 02 = 0.3 and k = 0. The associated queueing model 
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has the same parameters and 6% , but the other parameters are multiplied 
by n. The plots are shown without dividing by n. 

We ran the algorithm and the simulation for 50 time units. We used an 
Euler step of size h = 0.01, so we performed 5000 Euler iterations. In each 
Euler iteration we performed several iterations to calculate the matrix G in 
(11.1), which is used to calculate the instantaneous steady-state probability 

7Tl,2- 

Figures 1-4 show qi(t)/q2(t), ni^xit), q\{t) and z\${t) as functions of 
time t for a system initialized empty. After a short period, the pools fill up. 
Then qi(t) starts to grow, and immediately then fluid (customers) starts 
flowing to pool 2, causing z\^{t) to grow. Figures 1-4 show that, for practical 
purposes, steady state is achieved for t € [10,20]. 

In Figure 1 we see that once § 6 is hit, the ratio between the queues is kept 
at the target ratio 0.8. This is an evidence for the validity of the numerical 
solution, and a strong demonstration of the AP. In Figure 2 we see that 
initially, while qi = 0, 711,2 = 0. This lasts until £2,2 (i) + £1,2 (t) = m.2, at 
which time the space § is hit, specifically S b ), and the averaging begins. Once 
§ b is hit, 7Ti 2 becomes almost constant, even before the system reaches steady 
state. Thus the functions q\, q2 and z\ t 2 have exponential form, supporting 
the results of §9. 

We got x(t n ) = {qi(t n ),q 2 (t n ),z h 2{tn)) = (0.3639,0.4550,0.2385) and 
7Ti,2(ira) = 0.2 when the algorithm terminated. From (8.2), x* = (qf, q^, z\ 2 ) = 
(0.3667,0.4595,0.2375). From (8.3), we get tt* 12 = 0.2. 



Fig 1. ratio between the queues. Fig 2. 71-1,2 calculated at each iteration. 

Before solving the ODE, we can apply Theorem 10.2 to conclude that the 
solution will remain in A after it first hits A. 

12. Proof of Theorem 7.1. We have previously observed that the 
first two conclusions involving S + and S~ are valid. We now prove the three 
conclusions involving A, A + and A~. We will use the fact that a function 
mapping a convex compact subset of M m to W n is Lipschitz on that domain 
if it has a bounded derivative. Since we can always work with balls in M. m 
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Fig 3. trajectory of qi together with a 
simulated sample path of the stochas- 



FlG 4. trajectory of Zi t 2 together with 
a simulated sample path of the stochas- 



tic process Qi in a system initializing 
empty. 



tic process Z\^, in a system initializing 
empty. 



(which are convex with compact closure) , that in turn implies that a function 
mapping an open subset of M m to M n is locally Lipschitz whenever it has a 
bounded derivative on each ball in the domain; e.g., see Lemma 3.2 of [9]. 
The three sets A, A + and A - are convex. The key is what happens in A. 

For understanding, it is helpful to first verify this theorem in the special 
case r = 1, where the QBD process reduces to a BD process. Thus we first 
give a proof for that special case. 

Proof for the special case r = 1 . We use the fact that the ODE remains 
within A if it starts in A, so we are regarding A as an open connected convex 
subset of M 2 . The key component of the function ^ in A is 7Tx,2- We exploit 
the explicit representations in (6.18) and (6.21). From (6.1)-(6.4), the partial 
derivatives of A ± (x) and /J^(x) with respect to the three components of x, 
i.e., q±, q2 and zi^, are constants. From (6.18) and (6.21), we see that the 
partial derivatives of 71^2 (x) with respect to each of the three components 
of x exist, are finite and continuous. That takes care of A. 

We next consider A~ and A + ; the reasoning for these two cases is es- 
sentially the same, with (6.18) making it quite elementary. We see that 
7r i,2(x) — > and these partial derivatives approach finite limits as x — > 
Xb £ A - for x £ A, while 7Ti^(x) — > 1 and these partial derivatives ap- 
proach finite limits as x — > Xb £ A + for x £ A. In both cases we have a 
conventional heavy-traffic limit: p^{x) f 1 Hence, the partial 

derivatives of 711,2 (x) are continuous and bounded on § 6 . As a consequence, 
for any e-ball in S — §~ about x in A + , there exists a constant K such that 
1^1,2(^1) — 7Ti2(x2)| < -K"||xi — X2II3 for all x\ and X2 in the e-ball, where 
|| • H3 is the maximum norm onl 3 . A similar statement applies to A~. 

Hence we have completed the proof for r = 1. In closing, note that we 
cannot conclude that tti^{x) is even continuous on all of S, because for x G A 
we may have a sequence {x n : n > 1} with x n £ § + for all n (or x n £ S~ 
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for all n), with x n — > x as n — > oo, Kipfan) = 1 f° r & H n ( or = 0), while 
< 7Ti )2 (x) < 1. 

We now treat the general case. 

Proof of Theorem 7.1 in the general case. We first consider A. As in the 
case r = 1, we use the fact that the ODE remains within A if it starts in 
A, so we are regarding A as an open connected convex subset of R 2 . We 
will look at 7Ti2) an d thus the QBD, as a function of the variable x G A, 
which is an element of M 3 . By the definition of the matrices Aq, A\ and A2 
in (6.6) (see also the example in §6.2), these matrices are twice differentiable 
with respect to any of their elements. By the definition of the rates in (6.1)- 
(6.4), which are the elements of the matrices Aq, A\ and A2, these matrix 
elements in turn have constant partial derivatives with respect to each of 
the three real components of x at each x G A, i.e., with respect to qi, qi 
and z\ 2- It follows from Theorem 2.3 in He [7] that the rate matrix R in 
(6.15), which is the minimal nonnegative solution to the quadratic matrix 
equation Aq + RA\ + it! 2 ^2 = 0, is also twice differentiable with respect to 
the matrix elements of Aq, A\ and A2, and thus also with respect to the 
three real components of x at each x G A. 

It thus suffices to look at the derivatives with respect to one of the ele- 
ments of the matrices Aq, A\ or A2. It follows from the normalizing expres- 
sion in (6.16) and the differentiability of R, that ao is also differentiable. 
Hence, from (6.17), we see that wi^ is differentiable at each x G A, with 

(12.1) n' lt 2 = a' (I - R)~ 1 1+ + a (I - i?) _1 i?'(/ - J R) _1 1+. 
By differentiating (6.16), we have 

(12.2) a' (I - i?) _1 l + oq(I - R)- X B!{I - R^l = 0, 

so that a' is continuous. The continuity of R' and a' Q with respect to one of 
the elements of the matrices Aq, A± or A2 implies that the derivative Tr[ 2 
with respect to one of the elements of the matrices Aq, A± or A2 is finite 
and continuous on A, which in turn implies that the partial derivatives 
with respect to the three real components of x at each x G A are finite 
and continuous as well. Hence, VP is locally Lipschitz continuous on A, as 
claimed. 

We next show that 7ri 2 and thus ^ are locally Lipschitz continuous in 
neighborhoods of points in A + within S — S~ and of points in A~ within 
S — S + . We will only consider A + , because the two cases are essentially the 
same. In both cases, the situation is complicated starting from (12.1) because 
the entries of ao(x) become negligible, while the entries of (I — i?) _1 (x) 
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explode However, the two different limits cancel their effect. 

We exploit (6.19). The representation in (6.19) is convenient because now 
ao(x) —> oto(xb) as x — > Xb, where ao(xfc) is finite. All key asymptotics take 
place in R + . 

Since the crucial asymptotics involves only IR + , we see that we only need 
carefully consider one of the two regions, in this case the upper one. To 
obtain results about M + , from a process perspective, it suffices to replace 
the given QBD by a new QBD with the upper region and reflection at 
the lower boundary. The new QBD model involving only M + is equivalent 
to a relatively simple single-server queue. The net input is a linear com- 
bination of four Poisson processes, and so has stationary and independent 
increments. The queue length process in the revised model is an elementary 
MAP '/MSP/1 queue, as in §4 of [1], which has as QBD representation with 
rate matrix R + . 

For the asymptotics, the key quantities are the spectral radii of the matri- 
ces R + (x) and R~(x), say rj + (x) and r]~(x), and the way that these depend 
on the drifts 5+(x) and 5-(x) as x — > Xf,. The spectral radius r/ + (x) is the 
unique root in the interval (0,1) of the equation cfet[Aj~(x) + Af(x)r) + 
A^x) 1 )] 2 ] = 0, and similarly for rf~ (x); see (39) on p. 241 of [12], the 
Appendix of [13] and §4 of [1]. We see that r] + (x) — > rj + (xb) = 1 and 
r]~(x) — > 7]~(xb) < 1 as x — > x\, G A + . In general, we can represent powers 
of the matrix R (and similarly for R + and R~) asymptotically as 

(12.3) R n = vurf + o(r] n ) as n -> oo, 

where u and v are the left and right eigenvectors of the eigenvalue rj, re- 
spectively, normalized so that ul = 1 and uv = 1. Moreover, as r\ — > 1, the 
matrix inverse (/ — i?) _1 is dominated by these terms. 

Hence, we can do a heavy-traffic expansion of r/ + (x) and the related quan- 
tities as x — > Xb G A + with x G A, as in [2]; see the Appendix of [13]. As 
x — > Xb, all quantities in (6.19) have finite continuous limits as x — > Xb G A + 
except (I — R^lx))" 1 . We first have |<5+(a;)| — > and 8-(x) — > 6-(x(,), where 
< 5^{xb) < oo. We then obtain 

l- v + (x) = 4.z b )\S+(x)\+o(\6 + (x)\) 
(12 4) (7 - *(«)+)- = " + l X ^ b) +»((!- 1 + W)-) 

V + (xb)u+(x b ) . 

= c(x b )\6 + (x)\ + ° { ^ ) 

as x — >■ Xb and |(5+(x)| — > 0, where c, v + and u + are continuous functions 
of Xb on A + . The asymptotic relations in (12.4) together with (6.19) imply 
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that 

(12.5) |7ri )2 (x) - 7ri )2 (x6)| = ki,2(a;) - 1| = I - r{x)/{\ + r(x)\, 
where 

(12-6) r(x) ee a l\{~ fXl) ~ h(x b )\5 + (x)\ 

and |<5+(x)| — > 0, where h is a continuous function on A + . Hence, 
there exist constants K\ and AT 2 such that 

(12.7) Ki,2(z) - 7!"i,2(^b)| < #i|£+0e)| < AT 2 ||x - x b || 3 

for all x sufficiently close to Xb- Finally, we can apply the triangle inequality 
with (12.7) to obtain |7Ti j2 (xi) — 7ri i2 (x 2 )| < 2AT 2 ||xi — X2II3 for xi,x 2 in an 
e ball about x b in S — S - . Hence, 7Ti j2 (x) and thus ^ are locally Lipschitz 
continuous on A + within S — §~. Hence the proof is complete. □ 
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Appendix 

APPENDIX A: OVERVIEW 

In this appendix we present some supplementary material. In §B we an- 
alyze the system with an initial underloaded state. In that case we show 
that the approximating fluid models lead to our main ODE in finite time. 
In §C we elaborate on the algorithm in §11 and give two more numerical 
examples, including one where the solution, starting empty, first enters § 
in S + , and then moves from S + to A and then §>~~, with 7ri j2 experiencing 
a discontinuity. In §D we give some omitted proofs. Finally, in §E we draw 
conclusions and mention remaining open problems. 

APPENDIX B: TRANSIENT BEHAVIOR BEFORE HITTING § 

Recall that the FQR-T control is designed to respond to unexpected over- 
loads. We assume that the two classes operate independently until a time 
at which the arrival rates change, and the system becomes overloaded. Let 
be the time that the arrival rates change. We thus think of a system in 
steady state at time when the arrival rates change, with 

(B.l) 9 i(0) = 92 (0) = *i j2 (0) = *2,i(0) = °- 

In particular, qi(0) < k, and no sharing is taking place. A well-operated 
system tends to have a critically loaded fluid limit, yielding steady-state 
values 24,1(0) = mi and 22,2(0) = m,2, but we could also have an underloaded 
steady state, with 21,1(0) < mi and/or 22.2(0) < rri2 as well. 

The ODE in (4.1)-(4.2) can be regarded as the fluid limit of a sequence of 
overloaded queueing models. Class 1 was assumed to be overloaded due to 
the arrival rate being larger than the total service rate of service pool 1, while 
class 2 was overloaded either because its arrival rate was also too large (but 
less so than class 1), or because pool 2 was helping class- 1 customers. For the 
ODE, the system overload assumption translates into having 24,1 (i) = mi 
and 24,2 (£) + 22,2^) = m 2 for all t, so that the state space for the fluid 
limit was taken to be S. (The space § was defined in §5 and §7.) However, 
if either 21,1 (0) < m\ or 22,2(0) < m 2 , then the initial state is not in S, so 
we cannot use the ODE (4.1) to describe the system. There is a transient 
period [0,tg) during which the two service pools fill up, but the system is 
not yet overloaded. 

If sharing is eventually going to take place (i.e., if x* is in either A or § + ), 
then with initial conditions as in (B.l), we should certainly hit E> b . Sharing 
will begin only at a time T such that qi(T) — rg 2 (T) = k. In this section we 
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show that, if indeed x* 6 A U § + , then T < oo, where 



(B.2) 



T = w£{t > : x(t) G S b }. 



The transient period of the fluid system can be divided into two distinct 
periods: The first transient period, on the interval [0, T), lasts until the fluid 
limit hits S> b . The second transient period is the one starting at the hitting 
time T, and is described by the ODE (4.2). This period was analyzed in the 
previous sections. The first transient period is described by different ODE's, 
depending on the state of the system. These ODE's, for the initial condition 
in (B.l), are given in the proof of Theorem B.l below. 

We shall prove that T < oo under the extra assumption that at no time 
during [0, T) is 22,1 > 0. The assumption can be verified directly by solving 
the fluid model of the first transient period. We discuss this condition after 
the proof of Theorem B.l. 

Theorem B.l. If x* e AU §+, if (B.l) holds and if z 2 ,i(t) = for all 
t > 0, then T < oo, for T in (B.2). 

Proof. We start by developing the ODE to describe the system before 
hitting S. As before, we do not consider the original queueing model and 
prove convergence to the appropriate fluid limit, but instead we develop the 
ODE directly. We first consider the case in s 2 > (so that gf = 0), i.e., 
class 2 experiences no overload by itself (before pool 2 starts serving class- 1 
fluid). First, there is an initial period in which the pools are being filled with 
fluid. It is easy to see that as long as neither pool is full, the pool-content 
functions Za(t) behave as the fluid approximations for the number in system 
at time t in an M/M /oo queueing model with arrival rate Aj and service rate 
fj,i t i, i = 1, 2; e.g., see [14] (where it assumed that A = fx, so that A/// = 1). 
Therefore, the system evolution is described by the pair of ODE's 



These ODE's describe the dynamics of the two classes until one of the pools 
is full, i.e., until the time 



zi,i(t) = Ai - ^i,izi,i(t), 21,1(0) 

Z2,2(t) = A 2 - ^2,2^2,2 (*), 22,2(0) 



Ci 

C2 



and the unique solution to each ODE is 




(B.3) 



t\ = min inf{t > : Zi^it) 



mi}. 
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Since we assume that s 2 > 0, t\ is the time at which Z\ t \{b) = mi, and at 
this time we need to start considering q±. Clearly, q\ evolves independently 
of class 2 until qi(t) = k (when sharing is initialized). Let 

(B.4) t 2 = inf{t > ti : qi(t) = k}. 

Recall that k may be equal to 0, in which case t\ = t 2 . If t 2 > t±, then qi(t), 
t G [t±,t 2 ), evolves as the fluid approximation for the queue- length process 
in an Erlang-A model operating in the ED MS-HT regime, as in [21]. The 
ODE describing the evolution of q\ is 

(B.5) qi(t) = Ai - Hi^mi - 0iqi(t), h<t<t 2 , with q 1 (t 1 ) = 0, 
and its unique solution is 

gl (t) = Al "^ imi (l-e-^)), H<t<t 2 . 

Now, since qi(t 2 ) = k and q 2 {t 2 ) = 0, class- 1 fluid starts flowing to service 
pool 2, so that z\ 2 starts increasing. There is a time £3 such that, for t G 
[*2, ^3) , Qi(t) — K 1 Q2(t) = and all the excess class-1 fluid, that is not lost 
due to abandonment, is flowing to pool 2. Hence, z\ )2 satisfies the ODE 

zi,2(t) = (Ai - fj,i,xmi - Oik) - fj,i, 2 zi, 2 (t), t 2 <t<t 3 , with zi, 2 (t 2 ) = 0, 

whose unique solution is 

z 1>2 (t) = Al " ^ imi ~ ° 1K (l - e-^A , t 2 <t<t 3 . 
Mi,2 v ' 

Hence, t% = inf{t > t 2 : zi 2 (t) + z 22 (t) = m 2 }, so that at time t% both 
service pools are full, with qiit^) = k, q 2 (ts) = and qi(ts) — rq 2 (ts) = k. 
It follows that £3 is the time at which the fluid model hits the space §\ and 
the first transient period is over, i.e., t% = T for T in (B.2). 

Now we consider the second case in which q% > 0. In this case there are 
different scenarios: In the first scenario, pool 2 can be filled before pool 1, 
so that t\ = inf{t > : 212,2 = m 2 }, for t\ in (B.3). In that case q 2 begins 
to increase at time t\, evolving according to the ODE of the overloaded 
Erlang-A model 

92 (*) = A 2 - ^2,2^2 - 2 q 2 (t). 

However, by the assumption of the theorem, we have ruled out the case 
in which qi(t) — r 2 ^q 2 (t) = K 2> \, so that no class-2 fluid will flow to pool 
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1. Hence, from the beginning (time 0), Z\^\ increases until time tf > t\ at 
which Zi t i = m\. Then q\ increases, satisfying (B.5) with qi(tf) = 0. By the 
assumption on x* , and following Corollary 8.2, there exists a time T < oo 
such that qi(T) — rq2{T) = k. This is because rq2{t) < rq% < qf — k for 
all t <T. On the other hand, it follows trivially from the solution to (B.5), 
that qf is the globally asymptotically stable point of (B.5). Hence, for every 
e > 0, there exists t e such that q\(t) > qf — e for all t > t e . (This is because, 
by the initial conditions, qi(t) < qf for all £). Thus, we can find e > such 
that 

(B.6) rq% < qf - e - k < qi(t) - k for all t > £ e . 

The second scenario of the second case has pool 1 filled first at time t±, so 
that qi starts increasing according to (B.5). If q\ reaches K before qi starts 
increasing, then we have the same behavior as when > 0. However, if 
at time ti in (B.4) q2 > 0, then the two queues will continue increasing 
independently until time T. Once again, (B.6) can be shown to hold, so that 
T < oo. □ 

We can easily calculate the exact value of x(T) and use it to calculate 
the QBD drift rates 5+(x(T)) and 6-(x(T)) to find whether the positive- 
recurrence condition (6.10) holds at T, so that x{T) £ A. 

Remark B.l. (sharing in the wrong direction) In Theorem B.l we as- 
sumed that we never have £2,1 > 0. The reason is that, if 2:2,1 ever does 
become positive, then the fluid x never hits the region S. To see that this 
is so, suppose that for some time t& sharing is initialized, with class-2 fluid 
flowing to service pool 1. Then Z21 is increasing until a time t$ at which 
91(^5) — r q2{t^,) = k, and the AP begins to operate. At that time, Z2 t i will 
start decreasing according to the ODE 

Z2,l(t) = -fJQ,lZ2,l(t), t > £5, 

whose unique solution is 

(B.7) z 2 ,i(£) = *2,i(t 5 )e-^ l{t - t5) , £>£ 5 - 

Hence Z21 remains strictly positive for all t > £5, and § is never hit. 

Of course, the fluid state should be approaching a state in § as £ increases. 
However, if there is such a limit point, then that limit point itself typically 
will not be a stationary point, because if x did start at that limit point, then 
it will have to continue to move toward the final stationary point x*. 
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More generally, the failure of 22,1 to actually reach in finite time has 
practical implications for the FQR-T control in the original queueing system. 
It suggests that it should be beneficial to relax the one-way sharing rule, 
by introducing lower positive thresholds for z\p, and z%,i- For example, if 
z 2,i(t) > at some time t > 0, and at the same time sharing should be 
done in the other direction (because of a new overload incident, with class 

1 being more overloaded and needing to get help), then we will allow pool 

2 to start helping class 1, provided that 22,1 is smaller than some threshold 
S2,i > 0. In that case, if £2,1 (i) > ^2,1, then 2^1 will cross the threshold S2,i 
in finite time, as can be seen from (B.7). It remains to examine the system 
performance in response to such more complex transient behavior. 

For the cases covered by Theorem B.l, the system evolution over the 
entire halfline [0, 00) is a continuous "soldering" of the different ODE's, but 
at the soldering points ij, the functions under consideration are typically not 
differentiable. Hence, there is no single ODE that captures the full dynamics 
of the system. To see why, consider the case in which sJ, > and k > 0. 
Then, for t < t±, qi(t) = and q\ = 0, but for t\ < t < t%, qi(t) evolves 
according to (B.5), which typically has a strictly positive derivative at t±. 
Thus the left and right derivatives at t\ are not equal. Similar arguments 
hold for all the other soldering points. 

APPENDIX C: THE ALGORITHM AND MORE EXAMPLES 

C.l. More on the Algorithm. Let {t m : m = 0, 1,2, ...,ra} be the 
Euler steps, with t m+ \ — t m = h. In our experiments we found h = 0.01 to 
be a good candidate for the step size since it is small enough to minimize 
numerical errors, while the number of iterations needed for the ODE to reach 
its stationary point, is just a few thousands. Hence the algorithm takes only 
a few seconds to terminate. 

Let D{t) = qi(t) — 7-52 (i), denote the weighted difference between the 
two fluid queues. The discretization of the ODE in the numerical algorithm 
means that if, at step k — 1, D(tk_i) ^ S b but is close to it, then D(tk) may 
miss the boundary, even though the (continuous) ODE is at the boundary at 
time tfc. For that reason, if k — h < D(tk) < K + h, then we force x(ifc) to be 
in § 6 , by taking D(tk) = k- Once we have D(tk) = k we decide whether to 
keep staying on the boundary for the next Euler step, by checking whether 
(6.10) holds. According to the relation between the QBD drift rates at time 
tfc, we decide whether we should apply the AP, in order to find ^1^(^)5 or 
rather set nipitk) to zero or one. 

At any step in the algorithm, we must also decide which ODE to use. 
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That depends on the state of the system at each time, as described in §B. 
If the fluid state is not in §, as in the initial period of the example in §11 
and the example below, then we use the appropriate fluid model, as given 
in the proof of Theorem B.l. 

C.2. An Example with x* £ S~*~. We now consider the same example 
as in §11.3, except now we increase the arrival rate for class 1 substantially, 
so that x* 6 S + . In particular, we let Ai = 3.0 instead of 1.3. Once again, 
the system is initialized empty. That means that the fluid solution in S is 
moving between the two regions S b and S + . In particular, the solution first 
hits § & , as was proved in Theorem B.l, but it stays there for a short amount 
of time, and then crosses to S + . 

As before, we show the results multiplied by n = 1000 in the figures below. 
We see how Z2 : 2 starts increasing up to the time T in which £1,2 (T 1 ) +£2,2 (T 1 ) = 
m2- At this time Z2^(T) starts decreasing, and is replaced by class-1 fluid. 
Since no class-2 fluid is flowing to either of the service pool, all the class-2 
fluid output is due to abandonment. We can also observe that ^2,2 eventually 
hits 0, even though 212,2 satisfies the equation (B.7). This is due to the 
numerical errors, as described in §B. 

In steady-state we have q\ = X 2 /0 2 = 900/0.3 = 3000 and = (Ai - 
^lMi.i ~~ m 2^i, 2)/^ = 4000, as in Corollary 8.1 (ii). 




Fig 5. 21,2 when Ai exceeds the system's Fig 6. 22,2 when Ai exceeds the system's 
capacity. capacity. 

C.3. An Example With x* £ S~ Moving Through §+ And § b . 

The purpose of this example is to illustrate more complex dynamics. We 
make class 2 more overloaded than class 1, i.e., qf < q%, but we make 
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10 20 30 40 

Time 

Fig 7. qi when Ai exceeds the system's 
capacity. 




10 20 30 40 

Time 

Fig 8. q2 when Ai exceeds the system's 
capacity. 



the rates faster for class 1. Specifically, we considered the following model 
parameters: Ai = 13.0, A2 = 1.5, ^1,1 = 10.0, = 0.8, ^2,2 = 1 #1 = 2, 
02 = 0.2, r = 0.8 and k = 0. Note that the arrival, service and abandonment 
rates are all substantially greater for class 1 than for class 2. Nevertheless, 
class 2 is more overloaded than class 1: qf = 1.5 < 2.5 = gf. For this 
example, x* = (1.5,2.5,0) £ S~. 

We applied our algorithm to this example, letting the system start empty, 
i.e., x(0) = 0. The results are shown in the remaining figures, where here the 
results are scale up by multiplying by n = 100. Since the class- 1 arrival rate 
is so large, gi starts filling up rapidly, and becomes full first; see Figures 9 
and 10. Since k = 0, pool 2 starts helping class 1 as soon as pool 1 becomes 
full. At first, pool 2 has spare capacity. However, soon the spare capacity 
in pool 2 is exhausted. At that time, the solution hits §. Even at the time 
pool 2 becomes full, we have q\ > rq2, so that the solution enters S via S + , 
Thus pool 2 continues to help class 1 even after it is fully occupied, causing 
a dip in 2:2,2; see Figure 10. However, the ratio of the queue lengths qi/q2 
decreases from its peak of about 1.4 until it reaches the target ratio r = 0.8, 
producing the desired relation qi = rq 2 + k; see Figure 12 At that time 
(about t = 1.15, the solution that was in S + hits the set A. At that time, 
^i,2(x) jumps from 1 down to a value about equal to 0.6; see Figure 11. 
For an interval of time, the solution remains in A with the queue ratio fixed 
at the target r = 0.8. However, the load imbalance cause the solution to 
move within A, causing ^1,2 (x) to decrease until it reaches in the set A~, 
at about time t = 2.5. From A - , the solution moves immediately into S~, 
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where it rapidly converges to its stationary point. Of course, the stationary 
point x* is not actually reached in finite time. Indeed, after S~ is reached, 
z\ t 2 decreases exponentially to 0, but zi^it) > = z* for all t, consistent 
with Remark B.l. 

From the figures, it is evident that the numerical solution is not identical 
to the real solution. Because of the discrete step sizes in the Euler steps, the 
numerical solution misses A initially. In fact, we have to design the algorithm 
such that it "discovers" when § b is missed, and then force it to hit That 
is easy to do since, if x(t k ) £ § + and x(t k +i) £ S~, where t k is the time 
of the k th Euler step, k > 1, then §> b must have been missed. We can then 
compute qi(t k+1 ) and take rq 2 (t k+1 ) = qi(t k+1 ) - n. 

This discreteness of the numerical solution explains the erratic behavior 
of 7Ti 2 at the hitting time of A, shown in Figure 11. The thick vertical line 
just after time 1, exactly when r = 0.8 for the first time as can be seen from 
Figure 12, appears because -k\^ jumps between and 1 at each Euler step. 
These jumps are caused the solution missing S fc at first. 




10 20 30 40 10 20 30 40 50 



Fig 9. 21,2 when x'eS . Fig 10. 2:2,2 when x* G § . 

APPENDIX D: MISSING PROOFS 

PROOF of Theorem 7.2. Since < zi^ < m 2 and > in S, we only 
need to prove the upper bounds (7.5). For i = 1,2, let Ui(t) be the function 
describing the queue-length process (of queue i) in a modified system with 
no service processes (so that all the fluid output is due to abandonment). 
The queue-length process in the modified system evolves according to the 
ODE 

«i(t) = A; - d iUi (t), t > 0, 
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12 3 



Fig 11. 7ri,2 when x* G § over short Fig 12. ratio between the queues when 
initial interval. x* € S~ over short initial interval. 
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30 40 



40 50 



Fig 13. qi when x* 6 § 



FlG 14. q2 when x* € 



whose solution is 



Ui(t) 



+ «i(0) 



A, 



e~" l£ ,t > 0. 



It follows that itj(t) < Ui(0) V Aj/0j and, when Uj(0) = gi(0), the the right- 
hand side in (7.5) is an upper bound for Ui(t). We now show that this 
is also a bound for qi(t). For that purpose, define the auxiliary function 
fi(t) = qi(t) — Ui(t), t > 0, and observe that /i(0) = and fi{0) < 0. Hence, 
/ is decreasing at with f(t) < /(0) for all t G [0,(5) for some 5 > 0. This 
implies that qi(t) < U{(t) for all t 6 [0,5). 
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We now want to show that qi(t) < Ui(t) for all t > 0. For a proof by 
contradiction, assume that there exists some to > such that %(2o) > Ui(to), 
and let 

ti = sup{2 < t : qi(t) = Ui(t)}, 2 2 = inf{2 > 2 : %(2) = Ui(t)}. 

By the contradictory assumption and the continuity of q and u, we have 
< t\ < to < 2 2 . (2 2 may be infinite.) Then 

(D.l) qi(t) > Ui(t) for all t x < 2 < 2 2 . 

It follows from the mean- value theorem that there exists some £3 G (ii,io) 
such that 

m = *HW = ML > 0. 

2o — ti to — ti 
Hence, ^(23) > iiife). For i = 1, this translates to 

Ai-^imi-TTi^Ocfo)) [21,2(23)^1,2 + 22,2(23)^2,2] -0i<?i (23) > Ai-6»im(t 3 ). 
Thus, 

Qi (<?i (h) ~ ui(t 3 )) < -yui.imi - 7Ti >2 (x(2 3 )) [21,2(23)^1,2 + 22,2(23)^2,2] < 0, 

so that gi(2 3 ) < ui(2 3 ), contradicting (D.l). A similar argument holds for 
92- □ 

Proof of Corollary 8.1. If x* G S b , then the solution to (8.4) will have 
< z < m 2 , where the exact value of x* is readily seen to be the one in {£). 
If x* G § + , then q* — rq 2 > K, so that tt\ 2 = 1- Plugging 7r* 2 = 1 in the 
ODE for 21,2(2) in (4.2), we get 21,2(2) = 22,2(2)^2,2- Since at stationarity 
21,2(2) = 0, it follows that z 22 = 0, which implies that z\ 2 = m 2- Plugging 
this value of z* 2 , together with -k\ 2 = 1 when = 0, i = 1, 2, we get the 
values of and q 2 as in (m). 

Finally, if x* G §~, i.e., if — < k, then 7r£ 2 = 0, so that, by plugging 
this value of i:\ 2 i n the ODE for 2i, 2 (2) in (4.2), we see that ii j2 (2) = 
^1,221,2(2). Equating to zero, to get the value at stationarity, we see that 
z* 2 = 0. Plugging ^^2 = and z\ 2 = in the ODE for qi(t) and (? 2 (2), and 
equating these to zero, we get the values in (in). □ 

PROOF of Corollary 8.2. We prove (i) only. The proofs for (ii) and 
(Hi) are similar. First assume that x* G § b . Since 2 — 0> It follows from 
the expression for z* 2 in (i) of Corollary 8.1 that if q 2 > then qf — K> rq 2 . 
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If s 2 > then ql — n > fJ-i^s^/^i by Assumption A. For the other inequality 
we use the fact that 

* 0i0 2 (qf - k) - r0 1 (X 2 - [i 2 , 2 m 2 ) . 
*i,2 = a T~a m2 ' 

which implies the right-hand inequality in (8.8). 

Now Assume that (8.8) holds. It follows from the right-hand-side (RHS) 
inequality and the expression of z in (8.1) that 

= 0i02{qf - k) - r6i(\ 2 - H 2 , 2 m 2 ) 
r0ifi 2)2 + 2 A*i,2 
9 1 e 2 (r\ 2 /6 2 + m t2 m 2 /8i) - r6i{\ 2 - /x 2 2 m 2 ) 
- a V~q = m2 - 

From the left-hand inequality in (8.8), we see that, if s 2 = (and necessarily 
q 2 > = 4), then 

z > 0i02rq 2 - r0i(\ 2 - fi 2t2 m 2 ) _ 
r9ifJ, 2 , 2 + 6* 2 /xi,2 

If 4 > (and q% = 0), then 

z > 02Hl,2S 2 ~ r6 1 (X 2 - /i 2 ,2A 2 ) _ 02/^1,24 + r 6lV2,2S 2 _ ^ 
r9lfjL 2 ,2 + 6*2^1,2 rQ X ^2,2 + #2/«l,2 

Thus, if (8.8) holds, then s 2 < z < m 2 . This was shown to to imply that 
x* € in the proof of Theorem 8.1. (In fact, we have a stronger result, since 
we have z > s 2 . This is due to the requirement that qf — K,> /ii,2S 2 /0i, which 
is exactly Condition (/) in Assumption A.) 

We can show that the inequalities in (8.8) are strict if and only if x* € A 
by first observing that the inequalities are strict if and only if < z* < m 2 , 
and then directly calculate the QBD drift rates at the point x*. This is 
done in §10; see (10.5). It then follows that (6.10) holds at x* if and only if 
< z* < m 2 . □ 

APPENDIX E: CONCLUSIONS AND FURTHER RESEARCH 

In this paper we analyzed the deterministic ODE (4.1)-(4.2), arising as 
the MS-HT fluid limit of the overloaded X call-center model operating under 
the FQR-T control. In addition to being an interesting mathematical object 
in its own right, the ODE analyzed in this paper is a vital part of the 
FWLLN and FCLT in [17, 18]. We prove that the stationary point point 
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x*, which was developed heuristically in [15] using flow-balance arguments, 
is indeed the unique stationary point for the ODE. Moreover, we provided 
mild conditions under which the solution x(t) converges to x* as t — > oo. 
We also showed that the convergence to x* is exponentially fast, further 
justifying the steady-state analysis in [15]. 

We showed that the existence of a unique solution to the IVP (4.3) de- 
pends heavily on the characterization of the function ^ in (4.1) and its 
topological properties. These properties, in turn, depend on the state space 
of \&, and the regions of the state space in which \I' is continuous. These re- 
gions are further characterized by the probabilistic properties of the family 
of FTSP's {D t :t> 0}. 

To further relate to the model considered in our previous paper [15], 
in §B we considered the system at the time when the arrival rates first 
change. At that time, the system will typically be underloaded, so that the 
state space should not be S. After the change, we assume that the arrival 
rates are larger than the total service rate of the two pools. Specifically, we 
assumed Assumption A in §8. We then considered the first transient period 
[0, T), where T is the time at which S b is hit. Using alternative fluid models 
(ODE's), we showed that T < oo, under the conditions of Theorem B.l. 
The solutions to the fluid models during the first transient period are all 
exponential functions, so that this period also passes exponentially fast. 

Finally, we developed an efficient algorithm to solve the IVP (4.3), based 
on the matrix geometric method. This algorithm solves the different fluid 
models described in §B, and combines these solutions with the solution to 
(4.2) once the set A is hit, where the AP takes place. 

It remains to quantify or at least bound the number of times the fluid 
solution moves from one of the regions §+, §~ or S b to one of the others. 
Of course, the complexity of a solution is constrained by the fact that the 
solution path cannot cross over itself. It also remains to consider more com- 
plicated dynamics than provided by a single change in the arrival rates. 
The numerical algorithm applies more generally, but it remains to estab- 
lish mathematical results and examine the performance. For example, it 
remains to consider a second overload incident happening before the system 
has recovered from the first one. Finally, it remains to establish analogs of 
the results here for more complex models, e.g., with more than two classes 
and/or more than two service pools. 
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